/********************************
*
*                    FILE     FIT.C
*
*  This file contains functions necessary to fit a synthetic
*  light curve to an observed light curve.
*
************************************/

#include "header.h"


double Normalize( long nphases, double *phases, double *fluxes )
/*******************************
*
*  This function normalizes the synthetic light curve to preset values.
*
*  If mode.normalize is "SETVALUE": 
*      the synthetic light curve curve is normalized such that the flux 
*      is set equal to starsys.normflux at phase starsys.normphase.  Note
*      that the synthetic light curve might not have been calculated at 
*      starsys.normphase and, if not, the synthetic light curve is 
*      interpolated to starsys.normphase.
*
*  If mode.normalize is "MAXVALUE":
*      the synthetic light curve is normalized such that its maximum
*      value at any calculated phase is starsys.normflux.
*  
*  Input Data:
*        nphases   =  the number of data points in the light curve.
*        phases[]  =  the phases of the fluxes in array fluxes[].
*        fluxes[]  =  the fluxes in the synthetic light curve.
*
*  Output Data:
*        1)  The normalized light curve.
*        2)  The normalization factor.
*
**********************************/
{
   long int i;
   double interpflux, slope, maxflux, normfactor;

   if( nphases == 0 ) return ( 1.0 );

   interpflux  = -1.0;
   if( strcmp( mode.normalize, "SETVALUE") == 0 ) {
      for( i = 0; i < (nphases - 1); i++ ) {
         if( (starsys.normphase >= phases[i]  ) &&
             (starsys.normphase <  phases[i+1])   ) {
            slope = (fluxes[i+1] - fluxes[i]) / (phases[i+1] - phases[i]);
            interpflux = (starsys.normphase - phases[i]) * slope + fluxes[i];
         }
      }
      if( starsys.normphase == phases[nphases-1] )
         interpflux = fluxes[nphases-1];
      if( interpflux < 0.0 )
         quit("Normalization phase outside range of calculated phases.");
      if( interpflux == 0.0 )
         quit("Unable to normalize. Flux at normalization phase = 0");
      normfactor = starsys.normflux / interpflux;
   }
   else if( strcmp( mode.normalize, "MAXVALUE") == 0 ) {
      maxflux = fluxes[0];
      for( i = 1; i < nphases; i++) {
         if( fluxes[i] > maxflux ) 
            maxflux = fluxes[i];
      }
      if( maxflux == 0.0 )
         quit("Unable to normalize.  Maximum flux = 0.0");
      normfactor = starsys.normflux / maxflux;
   }
   else
     quit("Unrecognized normalization type in function Normalize().");

   for( i = 0; i < nphases; i++)
      fluxes[i] = fluxes[i] * normfactor;
   
   return( normfactor );
}


double lsqrNorm( long nphasesdata, double *phasedata, double *fluxdata, 
                        double *stdevdata, 
                 long nphasescalc, double *phasecalc, double *fluxcalc)
/*******************************
*
*  This function normalizes the synthetic light curve by minimizing
*  the variance of the fit of the synthetic light curve to the observed 
*  light curve.
*
*  The synthetic light curve will not generally be calculated at the same
*  phases as the observed light curve, so the synthetic light curve is 
*  interpolated at get its flux at the observed phases.
*
********************************/
{
   double slope, normfactor, weight, numerator, denominator, s[PHASEPOINTS];
   long int i, j;

   /***********************************
   *
   *  For every data point.......
   *
   *************************************/
   for( i = 0; i < nphasesdata; i++) {

      /****************************
      *
      *  First check that the phase of the data point falls 
      *  within the range of phases of the LSQR fit.
      *
      *****************************/
      if( (phasedata[i] >= mode.fitphase[0]) 
                               && (phasedata[i] <= mode.fitphase[1]) ) {

         /**********************
         *
         *  Then check that the phase of the data point falls within
         *  the range of phases covered by the synthetic light curve.
         *
         *************************/
         if( (phasedata[i] < phasecalc[0]) 
                            || (phasedata[i] > phasecalc[nphasescalc-1]) ) {
            printf("   Attempted LSQR fit at orbital phases beyond the ends");
            quit("of the synthetic light curve");
         }

         /**************************
         *
         *  Interpolate the synthetic light curve to the phase of
         *  the data point.
         *
         **************************/
         for( j = 1; j < nphasescalc; j++) {
            if( phasedata[i] == phasecalc[j-1] ) {
               s[i] = fluxcalc[j-1];
               j = nphasescalc;
            }
            else if( phasedata[i] == phasecalc[j] ) {
               s[i] = fluxcalc[j];
               j = nphasescalc;
            }
            else if( (phasedata[i] >= phasecalc[j-1])
                                        && (phasedata[i] < phasecalc[j]) ) {
               slope = ( fluxcalc[j] - fluxcalc[j-1] ) 
                                         /  (phasecalc[j] - phasecalc[j-1] );
               s[i] = slope * (phasedata[i] - phasecalc[j-1]) + fluxcalc[j-1];
               j = nphasescalc;
            }
            else 
               normfactor = 1.0;   /* do nothing */
         }
      }
   }

   numerator = 0.0;
   denominator = 0.0;
   for( i = 0; i < nphasesdata; i++) {
      if( (phasedata[i] >= mode.fitphase[0]) 
                             && (phasedata[i] <= mode.fitphase[1]) ) {
         weight = 1.0 / ( stdevdata[i] * stdevdata[i] );
         numerator = numerator + weight * fluxdata[i] * s[i];
         denominator = denominator + weight * s[i] * s[i];
      }
   }
   if( denominator <= 0.0 ) {
      printf("   Synthetic and observed light curves do not overlap");
      quit("in the LSQR phase range.  Unable to normalize.");
   }
   normfactor = numerator / denominator;

   for( i = 0; i < nphasescalc; i++) {
      fluxcalc[i] = fluxcalc[i] * normfactor;
   }

   return( normfactor );
}


double chiSqr( long nphasesdata, double *phasedata, double *fluxdata,
                  double *stdevdata, 
               long nphasescalc, double *phasecalc, double *fluxcalc)
/*******************************
*
*  This function finds the reduced Chi^2 of the fit of the synthetic light
*  curve to the observed light curve.  The Chi^2 is calculated only
*  over the phases within the range mode.fitphase[0] to [1] if
*  NORMALIZATION=  LSQR.  Otherwise it is calculated over the entire
*  light curve.
*
***********************************/
{
   double slope, x, chisquared, s[PHASEPOINTS], fitphase0, fitphase1;
   long int i, j, nchipoints;

   if( strcmp( mode.normalize, "LSQR")  == 0 ) {
      fitphase0 = mode.fitphase[0];
      fitphase1 = mode.fitphase[1];
   }
   else {
      fitphase0 = -0.5;
      fitphase1 =  1.0;
   }

   /***********************************
   *
   *  For every data point.......
   *
   *************************************/
   for( i = 0; i < nphasesdata; i++) {

      /****************************
      *
      *  First check that the phase of the data point falls 
      *  within the range of phases of the LSQR fit.
      *
      *****************************/
      if( (phasedata[i] >= fitphase0) 
                   && (phasedata[i] <= fitphase1) ) {

         /**********************
         *
         *  Then check that the phase of the data point falls within
         *  the range of phases covered by the synthetic light curve.
         *
         *************************/
         if( (phasedata[i] < phasecalc[0]) 
                   || (phasedata[i] > phasecalc[nphasescalc-1]) ) {
            printf("   Attempted to calculate Chi^2 at orbital phases");
            quit("beyond the ends of the synthetic light curve");
         }
         for( j = 1; j < nphasescalc; j++) {
            if( phasedata[i] == phasecalc[j-1] ) {
               s[i] = fluxcalc[j-1];
               j = nphasescalc;
            }
            else if( phasedata[i] == phasecalc[j] ) {
               s[i] = fluxcalc[j];
               j = nphasescalc;
            }
            else if( (phasedata[i] >= phasecalc[j-1])
                                  && (phasedata[i] < phasecalc[j]) ) {
               slope = ( fluxcalc[j] - fluxcalc[j-1] ) 
                                  /  (phasecalc[j] - phasecalc[j-1] );
               s[i] = slope * ( phasedata[i] - phasecalc[j-1] ) 
                                     + fluxcalc[j-1];
               j = nphasescalc;
            }
            else 
               x = 1.0;   /* do nothing */
         }
      }
   }

   chisquared = 0.0;
   nchipoints = 0;
   for( i = 0; i < nphasesdata; i++) {
      if( (phasedata[i] >= fitphase0) 
                   && (phasedata[i] <= fitphase1) ) {
         nchipoints++;
         x = ( s[i] - fluxdata[i] ) / stdevdata[i];
         chisquared = chisquared + x * x;
      }
   }
   /**************************
   *
   *  The choice of (nchipoints - 3) is a guess, based on fitting a
   *  minimum of 3 parameters: inclination, delta phi, and q.
   *
   *****************************/
   if( nchipoints > 3 ) {
      chisquared = chisquared / ( nchipoints - 3 );
   }
   else {
     chisquared = -1.0;
     printf("  Not enough data points to calculate Chi^2.");
   }

   return( chisquared );
}