
/**
 * FRETViewUtil.java
 * A simple class that provides utility methods for the FretView plugin
 * Created: June 15, 2006
 *
 * @author Nathan Stevens
 * @version 0.1 9/13/2006
 */
import java.util.*;
import java.text.*;

public class FRETViewUtil {
  public FRETViewUtil() { }
  
  // Method to return a DyeSpectra Object after passing a abs and emission spectra object
  public static DyeSpectra[] getDyeSpectra(SpectraData[] absa, SpectraData[] emsa) {
    int size = absa.length;
    DyeSpectra[] dsa = new DyeSpectra[size];
    
    // create the spectra data
    for(int i = 0; i < size; i++) {
      SpectraData abs = absa[i];
      SpectraData ems = emsa[i];
      
      // get the meta and header information
      String[] sa = abs.getMetaInfo();
      String[] dmeta1 = (String[])sa.clone();
      String[] dmeta2 = (String[])sa.clone();
      dmeta1[4] = "FRETView Data";
      dmeta2[4] = "FRETView Data";
      
      sa = abs.getHeader();
      String[] header1 = (String[])sa.clone();
      String[] header2 = (String[])sa.clone();
      
      // set the name and get the notes now
      String name = abs.getName();
      String notes = header1[5];
      
      // get the normalized data
      float[][] abs_data = getNormalizedData(abs.getData());
      float[][] ems_data = getNormalizedData(ems.getData());
      
      // create the abs and emission spectra data now
      SpectraData abs_sd = new SpectraData(name, "transmission", header1, dmeta1, abs_data, null);
      SpectraData ems_sd = new SpectraData(name, "Em", header2, dmeta2, ems_data, null);
      
      // create the new dye spectra object
      dsa[i] = new DyeSpectra(name, abs_sd, ems_sd, 1.0f, 1.0f);
      dsa[i].setNotes(notes);
    }
    
    return dsa;
  }
  
  // Method to return a normalize data array. 2D arrays only!
  public static float[][] getNormalizedData(float[][] xydata) {
    float[][] nxydata = new float[xydata.length][2];
    
    // find the maximum
    float max = 0.0f;
    for(int i = 0; i < xydata.length; i++) {
      if(xydata[i][1] > max) {
        max = xydata[i][1];
      }
    }
    
    // now normalize the data
    for(int i = 0; i < nxydata.length; i++) {
      nxydata[i][0] = xydata[i][0]; // x value
      nxydata[i][1] = (xydata[i][1]/max)*100.0f; // y value
    }
    
    return nxydata;
  }
  
  // Method to return xydata array in which a normalized absorption data as been
  // converted to
  public static float[][] getECData(float[][] xydata, float EC) {
    float[][] nxydata = new float[xydata.length][2];
    
    for(int i = 0; i < nxydata.length; i++) {
      nxydata[i][0] = xydata[i][0]; // x value
      nxydata[i][1] = (xydata[i][1]/100.0f)*EC; // y value
    }
    
    return nxydata;
  }
  
  // Method to convert to from nanometers to wavenumbers
  public static float nanometerToWavenumber(float nm) {
    float wavenumber = 1.0f/(nm*1.0e-7f);
    return wavenumber;
  }
  
  // Method to convert from nanometers to wavenumbers
  public static float[][] nanometerToWavenumber(float[][] xydata) {
    float[][] nxydata = new float[xydata.length][2];
    for(int i = 0; i < nxydata.length; i++) {
      nxydata[i][0] = 1.0f/(xydata[i][0]*1.0e-7f);
      nxydata[i][1] = xydata[i][1];
    }
    
    return nxydata;
  }
  
  // function to return the radiative lifetime given a dye spectra
  // this uses the Strickler-Berg formula. n is the refractive index
  // unfortunately this does work well.  It looks like I am missing a constant
  public static float calculateLifetime(DyeSpectra ds, float EC, float QY, float n, float[] range) {
    float lifetime = 0.0f;
    SpectraData abs_sd = ds.getAbsorption();
    SpectraData ems_sd = ds.getEmission();
    
    float[][] abs_data = getECData(abs_sd.getData(), EC);
    float[][] ems_data = ems_sd.getData();
    
    // set the range to integrate over now
    float minwl = abs_data[0][0] - 1;
    float maxwl = abs_data[abs_data.length - 1][0] + 1;
    
    if(range != null) {
      minwl = range[0];
      maxwl = range[1];
    }
    
    // set the dx value
    float dx = 1.0f;
    
    // calculate the first integral
    float int1 = 0.0f;
    for(int i = 0; i < (ems_data.length -1); i++) {
      float wl1 = nanometerToWavenumber(ems_data[i][0]);
      float wl2 = nanometerToWavenumber(ems_data[i + 1][0]);
      
      dx = wl1 - wl2;
      
      float x1 = ems_data[i][1]*wl1; // convert energy to relative quanta
      float x2 = ems_data[i + 1][1]*wl2;
      int1 += dx*(x1 + x2)/2.0f;
    }
    
    // calculate the second integral
    float int2 = 0.0f;
    for(int i = 0; i < (ems_data.length -1); i++) {
      float wl1 = nanometerToWavenumber(ems_data[i][0]);
      float wl2 = nanometerToWavenumber(ems_data[i + 1][0]);
      
      dx = wl1 - wl2;
      
      float v1 = (float)Math.pow(wl1, -3.0);
      float v2 = (float)Math.pow(wl2, -3.0);
      
      float x1 = (ems_data[i][1]*wl1)*v1;
      float x2 = (ems_data[i + 1][1]*wl2)*v2;
      int2 += dx*(x1 + x2)/2.0f;
    }
    
    // calculate the third integral
    float int3 = 0.0f;
    for(int i = 0; i < (abs_data.length -1); i++) {
      if(abs_data[i][0] >= minwl && abs_data[i][0] <= maxwl) {
        float v1 = nanometerToWavenumber(abs_data[i][0]);
        float v2 = nanometerToWavenumber(abs_data[i + 1][0]);
        
        dx = v1 - v2;
        
        float x1 = abs_data[i][1]/v1;
        float x2 = abs_data[i + 1][1]/v2;
        int3 += dx*(x1 + x2)/2.0f;
      }
    }
    
    // now calculate the kr
    float n2 = (float)Math.pow(n, 2.0);
    float kr = 2.88e-9f*n2*(int1/int2)*int3;
    
    // get the radiative lifetime now
    float t0 = 1.0f/kr;
    lifetime = QY*t0;
    
    // debug code
    System.out.println("DX " + dx);
    System.out.println("N2 " + n2);
    System.out.println("INT 1 " + int1);
    System.out.println("INT 2 " + int2);
    System.out.println("INT1/INT2 " + int1/int2);
    System.out.println("INT 3 " + int3);
    System.out.println("kr (s^-1)" + kr);
    System.out.println("t0 " + kr);
    System.out.println("Lifetime (s)" + lifetime);
    
    return lifetime;
  }
  
  // Method to calculate the forster distance
  public static float calculateR0(DyeSpectra donor, DyeSpectra acceptor, 
  float EC, float QY, float n, float k2) {
    float R0 = 0.0f;
    SpectraData abs_sd = acceptor.getAbsorption();
    SpectraData ems_sd = donor.getEmission();
    
    float[][] abs_data = getECData(abs_sd.getData(), EC);
    float[][] ems_data = ems_sd.getData();
    
    // find the index in the absorbtion data that corresponds to the first emission
    // wavelength
    int index = -1;
    int start = 0;
    for(int i = 0; i < ems_data.length; i++) {
      index = getWLIndex(abs_data, ems_data[i][0]);
      start = i;
      if(index != -1) {
        break;
      }
    }
    
    float dx = 1.0f;
    
    System.out.println("Index in absorption is " + index + "Wavelength " + ems_data[0][0]);
    
    // calculate the second integral
    float overlap = 0.0f;
    for(int i = start; i < (ems_data.length -1); i++) {
      float cm1 = ems_data[i][0]*1e-7f;
      float cm2 = ems_data[i + 1][0]*1e-7f;
      
      dx = cm2 - cm1;
      
      float x1 = ems_data[i][1]*cm1; // multiply by wavelenght to get the quanta
      float x2 = ems_data[i + 1][1]*cm2;
      overlap += dx*(x1 + x2)/2.0f;
    }
    
    // calculate the normalized overlap integral
    float norm_overlap = 0.0f;
    for(int i = start; i < (ems_data.length -1); i++) {
      if(index < abs_data.length - 1) {
        float cm1 = ems_data[i][0]*1e-7f;
        float cm2 = ems_data[i + 1][0]*1e-7f;
        
        dx = cm2 - cm1;
        
        float wl1 = (float)Math.pow(cm1, 4);
        float wl2 = (float)Math.pow(cm2, 4);
        
        float nFd1 = (ems_data[i][1]*cm1)/overlap; // normalized emission data
        float nFd2 = (ems_data[i + 1][1]*cm2)/overlap;
        
        float x1 = nFd1*abs_data[index][1]*wl1;
        float x2 = nFd2*abs_data[index + 1][1]*wl2;
        
        //System.out.println("I " + i + "Index " + index + " J = " + nFd1 + " K = " + x1);
        norm_overlap += dx*(x1 + x2)/2.0f;
        index++;
      }
    }
    
    // now calculate R0
    float n4 = (float)Math.pow(n, -4.0);
    float f1 = k2*n4*QY*norm_overlap;
    float exp = 1.0f/6.0f;
    float f2 = (float)Math.pow(f1, exp);
    R0 = 9.78e3f*f2;
    
    //* debug code
    System.out.println("DX " + dx);
    System.out.println("QY " + QY);
    System.out.println("N^-4 " + n4);
    System.out.println("K^2 " + k2);
    System.out.println("Overlap " + overlap);
    System.out.println("Normalized Overlap (M^-1 cm^3) " + norm_overlap);
    //System.out.println("f1 " + f1);
    //System.out.println("f2 " + f2);
    //System.out.println("Exp " + exp);
    System.out.println("R0 (Angs) " + R0); //*/
    
    return R0;
  }
  
  // function to get index of a wavelength given an array
  public static int getWLIndex(float[][] xydata, float wl) {
    int index = -1;
    
    float precision = (xydata[1][0] - xydata[0][0])/2.0f;
    for(int i = 0; i < xydata.length; i++) {
      float temp_wl = xydata[i][0];
      if(Math.abs(temp_wl - wl) <= precision) {
        index = i;
        break;
      }
    }
    
    //System.out.println("Wavelength " + wl + " found at index " + index + " Value " + xydata[index][0]);
    
    return index;
  }
  
  // function to return the multiplicative factor (MF) given a time scale
  public static float getTimeUnitMF(String unit) {
    float mf = 1.0f;
    if(unit.equals("picosecond")) {
      mf = 1.0e12f;
    }
    else if(unit.equals("nanosecond")) {
      mf = 1.0e9f;
    }
    else if(unit.equals("microsecond")) {
      mf = 1.0e6f;
    }
    else if(unit.equals("millisecond")) {
      mf = 1.0e3f;
    }
    else { // must be seconds
      mf = 1.0f;
    }
    
    return mf;
  }
  
  // function to return the time unit given a multiplicative factor (MF)
  public static String getTimeUnit(float mf) {
    String unit;
    if(mf >= 1.0e12f) {
      unit = "ps";
    }
    else if(mf >= 1.0e9f) {
      unit = "ns";
    }
    else if(mf >= 1.0e6f) {
      unit = "us";
    }
    else if(mf >= 1.0e3f) {
      unit = "ms";
    }
    else { // must be seconds
       unit = "s";
    }
    
    return unit;
  }
  
  // Method to calculate the energy transfer efficiency
  public static float[][] calculateEfficiency(float[] distances, float R0, float ratio) {
    // create the efficiency 2D matrix
    float[][] efficiency = new float[distances.length][2];
    
    // caculate the efficiency now
    float R06 = (float)Math.pow(R0, 6);
    for(int i = 0; i < distances.length; i++) {
      efficiency[i][0] = distances[i];
      
      float r6 = (float)Math.pow(distances[i], 6);
      efficiency[i][1]  = (ratio*R06)/(ratio*R06 + r6);
      //System.out.println("For \"r\"  = " + distances[i] + " E = " + efficiency[i][1]);
    }
    
    return efficiency;
  }
  
  // Method to calculate the energy transfer rate
  // kt  = E/(1-E)*kt
  public static float[][] calculateRate(float[][] efficiency, float tau) {
    float[][] transfer_rate = new float[efficiency.length][2];
    
    for(int i = 0; i < efficiency.length; i++) {
      transfer_rate[i][0] = efficiency[i][0];
      transfer_rate[i][1] = efficiency[i][1]/((1 - efficiency[i][1])*tau);
      //System.out.println("For \"r\"  = " + transfer_rate[i][0] + " kt = " + transfer_rate[i][1]);
    }
    
    return transfer_rate;
  }
  
  // Method to calculate fluorescence intensity of the donor
  // Fd = QYd*ECd(1-E)
  public static float[][] calculateDonorIntensity(DyeSpectra donorDye, float[][] efficiency, 
  float[] wavelengths) {
    float[][] donor_intensity = new float[wavelengths.length + 1][efficiency.length + 1];
    
    // get some variables
    SpectraData abs_sd = donorDye.getAbsorption();
    float[][] abs_data = getECData(abs_sd.getData(), donorDye.getEC());
    float QY = donorDye.getQY();
    
    // set the wavelength and distance axis
    for(int i = 0; i < wavelengths.length; i++) {
      donor_intensity[i+1][0] = wavelengths[i];
      for(int j = 0; j < efficiency.length; j++) {
        donor_intensity[0][j + 1] = efficiency[j][0];
      }
    }
    
    // calculate the fluorescence intensity now
    for(int i = 1; i < donor_intensity.length; i++) {
      float EC = getECAtWavelength(abs_data, donor_intensity[i][0]);
      
      for(int j = 1; j < donor_intensity[0].length; j++) {
        float Fd = QY*EC*(1 -  efficiency[j-1][1]);
        donor_intensity[i][j] = Fd;
        /*if(donor_intensity[i][0] >= 480 && efficiency[j-1][0] >= 150) {
          System.out.println("For Wl = " +  donor_intensity[i][0]  + " R = " + efficiency[j-1][0] + 
          " Fd " + Fd + " QY = " + QY + " EC = " + EC); // debug code
        }*/
      }
    }
    
    // now calculate the intensity value
    return donor_intensity;
  }
  
  // Method to calculate fluorescence intensity of the acceptor
  // Fa = QYa(ECa + E*ECd);
  public static float[][] calculateAcceptorIntensity(DyeSpectra donorDye, DyeSpectra acceptorDye, 
  float[][] efficiency, float[] wavelengths) {
    float[][] acceptor_intensity = new float[wavelengths.length + 1][efficiency.length + 1];
    
    // get some variables
    SpectraData abs_sd1 = donorDye.getAbsorption();
    float[][] abs_data1 = getECData(abs_sd1.getData(), donorDye.getEC());
    
    SpectraData abs_sd2 = acceptorDye.getAbsorption();
    float[][] abs_data2 = getECData(abs_sd2.getData(), acceptorDye.getEC());
    float QY = acceptorDye.getQY();
    
    // set the wavelength and distance axis
    for(int i = 0; i < wavelengths.length; i++) {
      acceptor_intensity[i+1][0] = wavelengths[i];
      for(int j = 0; j < efficiency.length; j++) {
        acceptor_intensity[0][j + 1] = efficiency[j][0];
      }
    }
    
    // calculate the fluorescence intensity now
    for(int i = 1; i < acceptor_intensity.length; i++) {
      float EC1 = getECAtWavelength(abs_data1, acceptor_intensity[i][0]);
      float EC2 = getECAtWavelength(abs_data2, acceptor_intensity[i][0]);
      
      for(int j = 1; j < acceptor_intensity[0].length; j++) {
        float Fa = QY*(EC2 +  EC1*efficiency[j-1][1]);
        acceptor_intensity[i][j] = Fa;
        //System.out.println("For Wl = " +  acceptor_intensity[i][0]  + " R = " + efficiency[j-1][0] + 
        //" Fa " + Fa); // debug code
      }
    }
    
    // now calculate the intensity value
    return acceptor_intensity;
  }
  
  // Method to calculate the total intensity
  public static float[][] calculateTotalIntensity(float[][] donor_intensity, 
  float[][] acceptor_intensity) {
    float[][] total_intensity = new float[donor_intensity.length][donor_intensity[0].length];
    
    float max = 0.0f; // used to normalized the data to 1
    for(int i = 1; i < total_intensity.length; i++) {
      total_intensity[i][0] = donor_intensity[i][0]; // set the wavelenght value
      for(int j = 1; j < total_intensity[0].length; j++) {
        if(i == 1) {
          total_intensity[0][j] = donor_intensity[0][j];
        }
        
        float total = donor_intensity[i][j] + acceptor_intensity[i][j];
        total_intensity[i][j] = total;
        
        // set the maximum
        if(total > max) {
          max = total;
        }
      }
    }
    
    // normalize the all data now
    for(int i = 1; i < total_intensity.length; i++) {
      for(int j = 1; j < total_intensity[0].length; j++) {
        donor_intensity[i][j] /= max;
        acceptor_intensity[i][j] /= max;
        total_intensity[i][j] /= max;
        
        //System.out.println("For Wl = " +  total_intensity[i][0]  + "r = " + total_intensity[0][j] + 
        //" Fd = " + donor_intensity[i][j] + " Fa = " + acceptor_intensity[i][j] + 
        //" Ft = " + total_intensity[i][j]); // debug code
      }
    }//*/
    
    return total_intensity;
  }
  
  // calculate the intensity
  public static float[][] calculateAcceptorLifetime(float[][] efficiency, float tau) {
    float[][] lifetime = new float[efficiency.length][2];
    
    for(int i = 0; i < lifetime.length; i++) {
      lifetime[i][0] = efficiency[i][1];
      lifetime[i][1] = (1.0f - efficiency[i][1])*tau;
    }
    
    return lifetime;
  }
  
  // function to return the EC for a particular wavelength
  public static float getECAtWavelength(float[][] xydata, float wl) {
    float EC = 0.0f;
    
    float precision = (xydata[1][0] - xydata[0][0])/2.0f;
    for(int i = 0; i < xydata.length; i++) {
      float temp_wl = xydata[i][0];
      if(Math.abs(temp_wl - wl) <= precision) {
        EC = xydata[i][1];
        break;
      }
    }
    
    return EC;
  }
  
  // Method to return the lifetime 
  public static float getLifetime(DyeSpectra dye) {
    float lifetime = dye.getLifetime();
    
    String unit = dye.getLifetimeUnit();
    float mf = getTimeUnitMF(unit);
    lifetime = lifetime/mf;
    
    return lifetime;
  }
  
  // function to format the lifetime based on the time
  public static String formatTau(float tau) {
    NumberFormat nf = NumberFormat.getInstance();
    DecimalFormat df = (DecimalFormat) nf;
    df.setDecimalSeparatorAlwaysShown(true);
    df.applyPattern("0.##");
    
    String taus = null;
    
    if(tau >= 1.0f) {
      taus = df.format(tau) + " s";
    }
    else if(tau >= 1e-3f) {
      taus = df.format(tau*1e3f) + " ms";
    }
    else if(tau >= 1e-6f) {
      taus = df.format(tau*1e6f) + " us";
    }
    else if(tau >= 1e-9f) {
      taus = df.format(tau*1e9f) + " ns";
    }
    else { // must be picoseconds
      taus = df.format(tau*1e12f) + " ps";
    }
    
    return taus;
  }
}//--end class FRETViewUtil
