/*****************************************************************************
#   Copyright (C) 1994-2008 by David Gordon.
#   All rights reserved.                           
#                                                                           
#   This software is part of a beta-test version of the Consed/Autofinish
#   package.  It should not be redistributed or
#   used for any commercial purpose, including commercially funded
#   sequencing, without written permission from the author and the
#   University of Washington.
#   
#   This software is provided ``AS IS'' and any express or implied
#   warranties, including, but not limited to, the implied warranties of
#   merchantability and fitness for a particular purpose, are disclaimed.
#   In no event shall the authors or the University of Washington be
#   liable for any direct, indirect, incidental, special, exemplary, or
#   consequential damages (including, but not limited to, procurement of
#   substitute goods or services; loss of use, data, or profits; or
#   business interruption) however caused and on any theory of liability,
#   whether in contract, strict liability, or tort (including negligence
#   or otherwise) arising in any way out of the use of this software, even
#   if advised of the possibility of such damage.
#
#   Building Consed from source is error prone and not simple which is
#   why I provide executables.  Due to time limitations I cannot
#   provide any assistance in building Consed.  Even if you do not
#   modify the source, you may introduce errors due to using a
#   different version of the compiler, a different version of motif,
#   different versions of other libraries than I used, etc.  For this
#   reason, if you discover Consed bugs, I can only offer help with
#   those bugs if you first reproduce those bugs with an executable
#   provided by me--not an executable you have built.
# 
#   Modifying Consed is also difficult.  Although Consed is modular,
#   some modules are used by many other modules.  Thus making a change
#   in one place can have unforeseen effects on many other features.
#   It may takes months for you to notice these other side-effects
#   which may not seen connected at all.  It is not feasable for me to
#   provide help with modifying Consed sources because of the
#   potentially huge amount of time involved.
#
#*****************************************************************************/
#include    <math.h>
#include    "dGetPrimerMeltingTemperature.h"
#include    "bool.h"
#include    "nNumberFromBase.h"
#include    "consedParameters.h"
#include    <iostream.h>



static int nEntropyOfDimer[4][4];


static int nEnthalpyOfDimer[4][4];                          




static bool bEntropyAndEnthalpyAlreadyInitialized = false;



static double dGetEntropy( char* szOligo, const int nLengthOfOligo ) {

   int nEntropy = 108;
   for( int nBase = 0; nBase < nLengthOfOligo - 1; ++nBase ) {
      int nFirstBase = nNumberFromBase[ szOligo[nBase] ];
      int nSecondBase = nNumberFromBase[ szOligo[ nBase + 1 ] ];
      nEntropy += nEntropyOfDimer[ nFirstBase ][ nSecondBase ];
   }

   return( - ((double)nEntropy) * 0.1 );
}



static double dGetEnthalpy( char* szOligo, const int nLengthOfOligo ) {

   int nEnthalpy = 0;
   for( int nBase = 0; nBase < nLengthOfOligo - 1; ++nBase ) {
      int nFirstBase = nNumberFromBase[ szOligo[nBase] ];
      int nSecondBase = nNumberFromBase[ szOligo[ nBase + 1 ] ];
      nEnthalpy += nEnthalpyOfDimer[ nFirstBase ][ nSecondBase ];
   }

   return( -((double)nEnthalpy) * 0.1);
}

      

   
static void initializeEnthalpyAndEntrophy() {

   nEntropyOfDimer[nA][nA] = 240; 
   nEntropyOfDimer[nA][nG] = 208; 
                           
   nEntropyOfDimer[nC][nA] = 129; 
   nEntropyOfDimer[nC][nG] = 278; 
                           
   nEntropyOfDimer[nG][nA] = 135; 
   nEntropyOfDimer[nG][nG] = 266; 
                           
   nEntropyOfDimer[nT][nA] = 169; 
   nEntropyOfDimer[nT][nG] = 129; 
                           
   nEntropyOfDimer[nA][nC] = 173;
   nEntropyOfDimer[nA][nT] = 239;
                          
   nEntropyOfDimer[nC][nC] = 266;
   nEntropyOfDimer[nC][nT] = 208;
                          
   nEntropyOfDimer[nG][nC] = 267;
   nEntropyOfDimer[nG][nT] = 173;
                          
   nEntropyOfDimer[nT][nC] = 135;
   nEntropyOfDimer[nT][nT] = 240;





   nEnthalpyOfDimer[nA][nA] = 91; 
   nEnthalpyOfDimer[nA][nG] = 78; 
                           
   nEnthalpyOfDimer[nC][nA] = 58; 
   nEnthalpyOfDimer[nC][nG] = 119;
                           
   nEnthalpyOfDimer[nG][nA] = 56; 
   nEnthalpyOfDimer[nG][nG] = 110;
                           
   nEnthalpyOfDimer[nT][nA] = 60; 
   nEnthalpyOfDimer[nT][nG] = 58; 


   nEnthalpyOfDimer[nA][nC] = 65;  
   nEnthalpyOfDimer[nA][nT] = 86;  
                            
   nEnthalpyOfDimer[nC][nC] = 110; 
   nEnthalpyOfDimer[nC][nT] = 78; 
                            
   nEnthalpyOfDimer[nG][nC] = 111; 
   nEnthalpyOfDimer[nG][nT] = 65; 
                            
   nEnthalpyOfDimer[nT][nC] = 56;  
   nEnthalpyOfDimer[nT][nT] = 91;  

}



double dGetPrimerMeltingTemperature( char* szOligo, const int nLengthOfOligo ) {

   if ( ! bEntropyAndEnthalpyAlreadyInitialized ) {
      bEntropyAndEnthalpyAlreadyInitialized = true;
      initializeEnthalpyAndEntrophy();      
   }
   

   double dLogDNA = 1.987 * log( 
      consedParameters::pGetConsedParameters()->dPrimersDNAConcentrationNanomolar_ / 4000000000.0 
                                );

   double dLogSalt = 16.6 * log10( 
      (double) consedParameters::pGetConsedParameters()->dPrimersSaltConcentrationMillimolar_ / 1000.0 
                                  );

//    cout << " oligo for enthalpy, entropy is: " ;
//    for( int nBase = 0; nBase < nLengthOfOligo; ++nBase)
//       cout << szOligo[nBase];
//    cout << endl;

//    cout << "enthalpy is " << dGetEnthalpy( szOligo, nLengthOfOligo ) << endl;
//    cout << "entropy  is " << dGetEntropy(  szOligo, nLengthOfOligo ) << endl;
//    cout << "dLogSalt = " << dLogSalt << endl;
//    cout << "fPrimersSaltConcentrationMillimolar_ = " <<
//       consedParameters::pGetConsedParameters()->fPrimersSaltConcentrationMillimolar_ << endl;


   double dMeltingTemperature = 
      ( dGetEnthalpy( szOligo, nLengthOfOligo ) * 1000.0 ) /
      ( dGetEntropy(  szOligo, nLengthOfOligo ) + dLogDNA )
      - 273.15 
      + dLogSalt;

//    cout << "melting temp = " << dMeltingTemperature << endl;
   return( dMeltingTemperature );
}