//Libraries
#include <stdio.h>
#include <math.h>
#include "mpfit.h"	//Downloaded from http://www.physics.wisc.edu/~craigm/idl/cmpfit.html
//#include "idl_export.h"	


//Declarations
int model(int m, int n, double *p, double *deviates, double **derivs, void *private);

struct vars_struct{
  double *x;
  double *y;
  double *ey;
};

//Removal routine
int sofieNO8_natural(	double *ptr_signal, 
			double *ptr_altitude, 
			double* ptr_time,
			char* ptr_mode, 
			double* ptr_residual, 
			int* ptr_data_size, 
			double* ptr_resChis, 
			double* ptr_maxChis, 
			int* ptr_chisFlag, 
			double* ptr_Aparameters,
			double* ptr_SigmaParameters,
			double* ptr_maxUncertainty,
			double* ptr_minUncertainty,
			int* ptr_bumpFlag, 
			double* ptr_maxBump,
			double* ptr_weight,
			double* ptr_minFitz,
			double* ptr_maxFitz,
			double* ptr_startBA,
			double* ptr_endBA)
{

/***********************************************************************************
*
* This function removes the thermal oscillation from the signal data array (ptr_signal). 
* See the function "model" at the end of this file for details on the model definition
* used.
* The extinction is first calculated, taking the measurement at z0 = 140 km as the exoatmospheric reference. 
* It is assumed for this calculation that above 140km no attenuation due to the atmosphere is present. Therefore,
* above z0 the oscillation in the extinction is completely exoatmospheric. The oscillation is fitted with the
* model above 140km and extrapolated over the entire data. The residual (extinction data - fitted oscillation)
* is calculated and returned to the address indicated by "ptr_residual". The model's decay and frequency are
* fixed parameters in the model definition. 
*
* Required Input: (All inputs are pointers to arrays containing the specified data)
*
*   ptr_signal............original radiometer signal array (counts) 
*   ptr_altitude..........altitude for the event measurement (km) (convention followed is that lowest index is lowest altitude)
*   ptr_time..............time for the event measurement (sec) 
*   ptr_mode..............occultation mode, 0 = sunrise, 1 = sunset 
*   ptr_residual..........address where the residual oscillation is stored 
*   ptr_data_size.........number of datapoints in the event
*   ptr_resChis...........address where the resulting reduced Chi^2 parameter will be saved
*   ptr_maxChis...........threshold for the maximum expected reduced Chi^2 parameter (normally 3-5 are acceptable thresholds since this is an adhoc calculation)
*   ptr_chisFlag..........address where the flag indicating if the reduced Chi^2 parameter complied with the specified threshold will be stored
*   ptr_Aparameters.......address where the array for the resulting parameters will be stored
*   ptr_SigmaParameters...address where the array for the 1 sigma uncertainties of the parameters will be stored
*   ptr_maxUncertainty....address where the array for the maximum uncertainty at each altitude will be stored
*   ptr_minUncertainty....address where the array for the minimum uncertainty at each altitude will be stored
*   ptr_bumpFlag..........address where the flag indicating if there is negative extinction in the residual will be stored
*   ptr_maxBump...........threshold for the maximum expected negative extinction (in extinction units) (example: -0.0003, -0.0002)
*   ptr_minFitz...........altitude at which the fit starts (in km)
*   ptr_maxFitz...........altitude at which the fit stops (in km)
*   ptr_startBA...........time at which the Balance Adjustment (BA) procedure starts (in s). Set to 0.0 if the function should calculate the start of the
*		..........BA procedure. 
*   ptr_endBA.............time at which the BA procedure ends (in s). Set to 0.0 if the function should calculate the end of the BA procedure.
*		..........WARNING: The given BA position must be outside the BA, otherwise the fit will not converge (it will take data during the BA procedure).
*
* Optional input: 
*
*   ptr_weight............uncertainty in measurements for the fit (error of each measurement). If the input is 0, the default estimated value is used (0.504).
*
* Output:
*   ptr_residual..........contains the residual signal after applying correction (counts)
*   ptr_resChis...........contains the resulting reduced Chi^2 for the fit
*   ptr_chisFlag..........address where the flag indicating if the reduced Chi2 parameter complied with the specified threshold will be stored
*			  0 if complied, 1 if not
*   ptr_Aparameters.......address where the array for the resulting parameters will be stored
*   ptr_SigmaParameters...address where the array for the 1 sigma uncertainties of the parameters will be stored
*   ptr_maxUncertainty....address where the array for the maximum uncertainty at each altitude will be stored
*   ptr_minUncertainty....address where the array for the minimum uncertainty at each altitude will be stored
*   ptr_bumpFlag..........address where the flag indicating if there is negative extinction in the residual will be stored
*			  0 if complied, 1 if not.The bump is only searched at altitudes below the "beginning of BA", around 180km.

*** IMPORTANT REMARKS ***
	- Memory space allocation for each of the variables indicated by the pointers must be done prior to the call of the function. 
	- Also scalar values must have proper memory space allocation, since they are treated as pointers.
	- The reduced chi^2 parameter is an ad hoc calculation. The error in each data point was estimated. Appropiate weights can be provided as a parameter.
	- Time is assumed to increase with index for both sunset (SS) and sunrise (SR) events

* Authors: David Gomez Ramirez, Hampton University
*          Dr. John McNabb, Hampton University
*	   
************************************************************************************/
/*********CHANGE LOG************
*Version 2 - New model incorporated. The model can now take into account the data after the BA procedure.
	   - The new values for the decay and frequency introduced.  
*Version 3 - Performs fits accross the BA in the calibration events. Also works for real events with no BA. 
	   - The fitting range (max and minimum altitude in km) is now a parameter. 
*Version 4 - Implements an algorithm to identify data including pre and post BA procedure 
*Version 5 - Implements an invert function for automatic retrieval of BA in sunset or sunrise events
*Version 6 - Start and end of BA added as parameters.
*Version 7 - No inverting the data. Bump flag corrected, provided BAstart and BAend as parameters. 
*Version 8 - BA start and end can be provided or calculated independently. 
/*******************************/

	//	printf("Initializing function. \n");
//Constants
	const double d_pi = 3.14159265;

	const double i_z0 = 140.0;					//Altitude (km) at which the fitting starts
	const float f_zStep = 0.2;					//Altitude step from data point to data point (Km)	
	const float f_peakThreshold =  15.0;				//Change between two consecutive normalized points 
									//(in the finite difference calculation) above which it is considered a peak. 
	
	const int i_nParameters = 9;					//Number of parameters in the model
	const int BAoffset = 3;						//After BA init and BA end are found, the offset in datapoint to safely avoid the BA. 
	const int dataPointsPeakThreshold = 5;				//Number of data points after a peak is detected that are ignored for 
									//the next peak detection (to avoid the noise in the positive or negative edge of 
									//the BA counting as successive peaks). 

	//initialGuess = [amplitude, decay, frequency, phase angle, slope or linear term] Look at model for description.
	const double initialGuess[5] = {0.00004, -0.029877, 0.662752, 0.0, 0.0};	//Initial guess for the last 5 parameters, corresponds to the oscillation.
	int posGivenBAcorrection = 20;


//Variable declaration
	double darr_extinction[*ptr_data_size];
	double darr_spikeDetection[*ptr_data_size];
	double darr_timeDif[*ptr_data_size];
	double darr_extrapolated[*ptr_data_size];

	double d_minFitz = *ptr_minFitz;				//Min altitude for the fit (Km)
	double d_maxFitz = *ptr_maxFitz;				//Max altitude for the fit (Km)
	
	double darr_ey[*ptr_data_size];					//Weight = 1/variance = 1/dev^2
	//Initialize array to default value
	int i;
	for(i = 0; i <= (*ptr_data_size)-1; i++)
	{
		if(*ptr_weight == 0.0){
			darr_ey[i] =0.5407;				//Default estimated error in data. According to Gordley et al. in "The solar occultation 
									//for ice Experiment" article, doi> 10.1016/j.jastp.2008.07.012 the error is below the
									//digitalization level. This means that the uncertainty in the measurement must be around
									//half 1 count or 0.5 counts.   
		}else{
			darr_ey[i] =*ptr_weight;			//Use provided error for the data
		}
	}

	int b_calculateBAend = 0;					//Logic variable to identify whether to calculate the BA end or not
	if(*ptr_endBA == 0.0){
		b_calculateBAend = 1;
	}

	int b_calculateBAstart = 0;					//Logic variable to identify whether to calculate the BA start or not
	if(*ptr_startBA == 0.0){
		b_calculateBAstart = 1;
	}
	

	double d_x0;	
	double d_t0;
	int index_z0 = 0;
	int index_minFit;
	int index_maxFit;
	int index_BAInit;
 	int index_BAEnd;
 
//Variables required for the mpfit function
	struct vars_struct vars;						//Contains the data, x,y and ey (weight)	

	double perror[i_nParameters];						//Error array for the resulting parameters
	for(i = 0; i < i_nParameters; i++)
	{
		perror[i] = 0.0;						//Initialize the array 
	}


	mp_par pars[i_nParameters];						//Array containing struct with restrictions
	bzero(pars, sizeof(pars));           					//Initialize constraint structure

	int i_statusMPFIT;							//Stores the status of the fit

//Variables for the deviation calculation
	int iarr_deviationScale[] = {1,1,1,1,1,1,1,1,1};			//Scale factor for the sigma uncertainty of each parameter
	double maxFit[*ptr_data_size];
	double minFit[*ptr_data_size];

	double errorTau = 0.001673;						//Use the std. deviations as uncertainty for these parameters
	double errorOmega = 0.001947;						//These values were obtained when fitting the calibration events. 
										//The average decay and frequency is fixed for the fits of the 
										//real events. The standard deviation of those mean values are 
										//taken as the error of the parameter. 

//Variables for the type of event
	int mode = *ptr_mode;


//Variables to temporarily store index data.
	int temporalIndex = 0;


//////////////////////////////////////////////
///////////////PROGRAM START//////////////////
//////////////////////////////////////////////

//Extract index for altitude Z0 = 140km
	for(i = 0; i < (*ptr_data_size); i++)
	{
		if((float)ptr_altitude[i] >= i_z0 && (float)ptr_altitude[i] < (i_z0+f_zStep)){				//When altitude exceeds the exoatmospheric threshold
			index_z0 = i;
				//	printf("Altitude Z0 %f.\n", ptr_altitude[i]);
			break;		
		}
	}	
		//	printf("Obtained Z0 %i, %i.\n", index_z0, *ptr_data_size);

//Extract signal and time values at Z0 = 140km
	d_x0 = ptr_signal[index_z0];
	d_t0 = ptr_time[index_z0];
		//	printf("Calculated X0 and T0.\n");

//Calculate time differential (tdifferential = time - time at Z0)		
	for(i = 0; i < *ptr_data_size; i++ )
	{
		darr_timeDif[i] = (ptr_time[i]) - d_t0;    
	}	
		//	printf("Calculated time differential.\n");
		

//Get fitting index
	for(i = 0; i < *ptr_data_size; i++)
	{
		//This algorithm is used to extract the index at which the min and max fit altitude occurs. The f_zStep indicates the difference between one altitude and the next, so that the indexes are extracted at exactly the altitude required. This was implemented for the events in which the altitude is estimated.
		if((float)ptr_altitude[i] >= d_minFitz && (float)ptr_altitude[i] < d_minFitz+f_zStep){	//Altitude truncated to floating precision for comparison
			index_minFit = i;
				//	printf("Altitude Min Fit %f. Index Min Fit %i \n", ptr_altitude[i], i);
		}		
		if((float)ptr_altitude[i] >= d_maxFitz && (float)ptr_altitude[i] < d_maxFitz+f_zStep){	//Altitude truncated to floating precision for comparison
			index_maxFit = i;
				//	printf("Altitude Max Fit %f. Index Max Fit %i \n", ptr_altitude[i], i);
		}
	}	

		//	printf("Fitting index obtained.\n");

//In a SS event, altitude decreases with index, so index_minFit and index_maxFit need to be inverted.
//*ptr_mode: 0 = sunrise, 1 = sunset 
	if(*ptr_mode){
		temporalIndex = index_minFit;
		index_minFit = index_maxFit;
		index_maxFit = temporalIndex;
	}


//Extract data for the fit
	int i_fitElements = (index_maxFit-index_minFit) + 1;				//The +1 is to account for the element AT maxFit

	double darr_timeDifFit[i_fitElements];
	double darr_signalFit[i_fitElements];
	double darr_eyFit[i_fitElements];

	for(i=0; i<i_fitElements;i++)
	{
		darr_timeDifFit[i] = darr_timeDif[index_minFit+i];
		darr_signalFit[i] = ptr_signal[index_minFit+i];
		darr_eyFit[i] = darr_ey[index_minFit+i];				
	}

		//	printf("Fit variables initialized.\n");


//Balance Adjustment (BA) procedure detection.
	//The initial value for the BA end and init is the index of the maximum altitude (in order to get the sizes of the arrays correctly). 
	//This depends on the type of event
	//*ptr_mode: 0 = sunrise, 1 = sunset 
	if(*ptr_mode){	//SS
		index_BAEnd = 0;
		index_BAInit = 0;

	}else{		//SR
		index_BAEnd = i_fitElements-1;
		index_BAInit = i_fitElements-1;
	}

	int b_initBAfound = 0;
	int b_endBAfound = 0;
	int counterDatapoints = 0;			//Counts for a certain amount of datapoints after the beginning of the BA so that the end of BA is correctly identified

//Verify if the BA positions were provided
	//Look for the corresponding index, assuming that time always increases with index (which is the case for both SS and SR events) 
	for(i = 0; i<i_fitElements;i++){
		//Sunset or sunrise makes a difference in this case. *ptr_mode: 0 = sunrise, 1 = sunset 
		if(ptr_time[index_minFit + i] < *ptr_startBA && !*ptr_mode && !b_calculateBAstart){
			index_BAInit = i;
		}	
		if(ptr_time[index_minFit + i] < *ptr_endBA && !*ptr_mode && !b_calculateBAend){
			index_BAEnd = i;	
		}	
		if(ptr_time[index_minFit + i] < *ptr_startBA && *ptr_mode && !b_calculateBAstart){
			index_BAInit = i;		
		}	
		if(ptr_time[index_minFit + i] < *ptr_endBA && *ptr_mode && !b_calculateBAend){
			index_BAEnd = i;		
		}	
	}	//End for

	//Look for the BAend or BAinit that remains. 
	for(i = 0; i<=i_fitElements-2;i++){		//Do not take the last differential (because the resultant differential compares against zero which 
							//incorrectly looks like a peak).
 
		darr_spikeDetection[i] = darr_signalFit[i+1] - darr_signalFit[i];	//Finite difference

		if((fabs(darr_spikeDetection[i]) > f_peakThreshold) && !b_initBAfound){	
		//Absolute value is to avoid distinction between SS and SR events (spikes are inversed)
		//Notice however that index_BAInit in a SS is actually the end of the BA procedure following altitude. 
		
		//Verify if the BAinit was already found
			if(b_calculateBAstart){
				index_BAInit = i;	//Beginning of BA following index.
			}
			temporalIndex = i;		// We save this index for the special case of a SS event with fitmax during the BA procedure with BAstart given.
							// In this special case, the start of the BA adjustment is given, but it is not part of the fitting data (it occurs
							// after fitMax). When the program looks for the BA end (which we are actually interested in) it will confuse it with
							// BA init because it will only see one spike. Since BA init was given it will not save the value, so we 
							// use this temporal index. The index_BAinit is correctly set given conditions after this loop. Notice if BAend
							// is given, the first spike will look to the code as BAinit, when it is actually BAend in this special case.
							// That is also taken into account in the logic after the loop. 
			b_initBAfound = 1;	
			continue;
		}	

		//This counter prevents that successive peaks are taken into account. Only first peak must be considered.  
		if(b_initBAfound && counterDatapoints < dataPointsPeakThreshold+1){
			counterDatapoints++;
		}

		if((fabs(darr_spikeDetection[i]) > f_peakThreshold) && (counterDatapoints > dataPointsPeakThreshold)){

			if(b_calculateBAend){
				index_BAEnd = i;	//End of BA following index.
			}
			b_endBAfound = 1;
		}	
	}	//End for




//Apply BA offset to make sure not a single datapoint fo the BA is taken into account for fit. The offset depends on the type of event (SS, SR)
	if(b_endBAfound){	//Implies BAinit was also found. Here the logic is INDEPENDENT of the type of event	
		index_BAEnd += BAoffset;		
		index_BAInit -= BAoffset;
	}

	if(b_initBAfound && !b_endBAfound){	//The logic DEPENDS on the type of event
	//*ptr_mode: 0 = sunrise, 1 = sunset
		if(*ptr_mode){	//SS
			
			if(b_calculateBAend){	
				//Save the real position of BAend (this applies if BAinit was given).
				index_BAInit = temporalIndex + dataPointsPeakThreshold + BAoffset; 
			}else{
				index_BAInit = index_BAEnd + BAoffset;			//This applies if BAend was given. 
				index_BAEnd = 0;
			}

		}else{		//SR
			index_BAInit -= BAoffset;
		}
	}

	//Create new arrays with the data without the BA
	//The minus one is to correctly fill the array assuming it is counting from 0. 
	int fitElements_woBA =  i_fitElements - abs(index_BAEnd - index_BAInit)-1;	//Abs to avoid issues for the SS events which 
											//have a higher index for BAInit than BAEnd.
	double darr_signalWithoutBA[fitElements_woBA];
	double darr_timeWithoutBA[fitElements_woBA];
	double darr_altitudeWithoutBA[fitElements_woBA];
	
//Take out the BA region for the fit.
	//This will also depend on the type of the event and on whether an end and beginning of BA were found.  

	for(i = 0; i<=i_fitElements-1;i++){
	if(b_endBAfound || !*ptr_mode){	//This logic still holds for SR events but not for SS events
		//Take out the BA or the portion of it in the data
		if(	i < index_BAInit){
			darr_signalWithoutBA[i] = darr_signalFit[i];
			darr_timeWithoutBA[i] = darr_timeDifFit[i];
		}else{
		if(	i > index_BAEnd){
			darr_signalWithoutBA[i- (index_BAEnd - index_BAInit) - 1] = darr_signalFit[i];
			darr_timeWithoutBA[i- (index_BAEnd- index_BAInit) - 1] = darr_timeDifFit[i];
		}}
	} else {//No end BA found and the event is a SS. 
		//Take out the BA which is after the detected spike. 
		if(	i > index_BAInit){
			darr_signalWithoutBA[i-index_BAInit-1] = darr_signalFit[i];
			darr_timeWithoutBA[i-index_BAInit-1] = darr_timeDifFit[i];
		}
	}//end else
	}//end for
		//	printf("BA removal procedure completed.\n");


//Data and variables definition for the mpfit function
	vars.x = darr_timeWithoutBA;
	vars.y = darr_signalWithoutBA;
	vars.ey = darr_eyFit;

	mp_result result;
	bzero(&result,sizeof(result));		//Zero results structure, (initialization)
	result.xerror = perror;
	
	//Format of the parameters = [x0, Cba pre, Cba post, indexBA, Amp, Decay, Frequency, Phase, Slope]
	double parameterValues[] = {d_x0, 1, 1, index_BAInit, initialGuess[0], initialGuess[1], initialGuess[2], initialGuess[3], initialGuess[4]};	


	if(!b_endBAfound){					//No complete BA adjustment.
		parameterValues[3] = -1.0;			//Initial guess. These only apply to B16	
		//	printf("No end of BA found \n");
	}

	pars[0].fixed = 1;                      //Fix parameter x0
	pars[3].fixed = 1;                      //Fix parameter BA time
	pars[5].fixed = 1;                      //Fix parameter Decay
	pars[6].fixed = 1;                      //Fix parameter Frequency

	pars[1].fixed = 0;                    
	pars[2].fixed = 0;                     
	pars[4].fixed = 0;                    
	pars[7].fixed = 0;                      
	pars[8].fixed = 0;                     


//Configuration for MPFIT
	mp_config config;
	bzero(&config, sizeof(config));
    config.maxfev = 1000;
    config.nofinitecheck=1;

	//	printf("MPFIT initialized.\n");

//Calling fitting routine
	i_statusMPFIT = mpfit(model, fitElements_woBA, i_nParameters, parameterValues, pars, &config, (void *) &vars, &result);
	//	printf("Fit performed.\n");

//Retrieve Red Chi^2 for the fit:
	*ptr_resChis = result.bestnorm/(result.nfunc-result.nfree);	//Chi^2 / DoF is the reduced chi^2 parameter needed

//Retrieve final value for parameters
	for(i = 0; i<i_nParameters; i++)
	{
		ptr_Aparameters[i] = parameterValues[i];
		//	printf("Parameter %f",parameterValues[i]);
	}

//Extrapolation
	double poly,s,e, d_CbaPre, d_CbaPost, d_indBA, d_amp,d_tau,d_ome,d_phi,d_slo;
	d_x0 = parameterValues[0];
	d_CbaPre = parameterValues[1];
	d_CbaPost = parameterValues[2];
	d_indBA = parameterValues[3];
	d_amp = parameterValues[4];
	d_tau = parameterValues[5];
	d_ome = parameterValues[6];
	d_phi = parameterValues[7];
	d_slo = parameterValues[8];

	/*Step functions*/
	double darr_stepPre[*ptr_data_size];
	double darr_stepPost[*ptr_data_size];

	for(i = 0; i<*ptr_data_size ;i++){

		if((i < d_indBA+index_minFit)|| (d_indBA < 0)){ 		//If the indBA is negative, it means that there is no data post BA
			darr_stepPre[i] = 1;
			darr_stepPost[i] = 0;
		}else{
			darr_stepPre[i] = 0;
			darr_stepPost[i] = 1;
		}
	}

	//Evaluating model
	for (i=0; i<*ptr_data_size; i++) {
		poly=-d_amp*sin(d_phi)+d_slo*darr_timeDif[i];
		s=sin(darr_timeDif[i]*d_ome+d_phi);
		e=exp(darr_timeDif[i]*d_tau);
		darr_extrapolated[i] = (d_CbaPre*darr_stepPre[i] + d_CbaPost*darr_stepPost[i])*d_x0*(1 - (d_amp*e*s+poly));
		
		//Residual calculation
		ptr_residual[i] = ptr_signal[i] - darr_extrapolated[i] + d_x0;  
	}


	//	printf("Extrapolation performed.\n");


//Complete uncertainties for fixed parameters (Tau and Omega) with their calculated standard deviations
	result.xerror[5] = errorTau;
	result.xerror[6] = errorOmega;
	
//Return the 1 sigma uncertainties of the parameters
	for(i = 0; i<i_nParameters; i++)
	{
		ptr_SigmaParameters[i] = result.xerror[i];
	}
		

//Initialize maxFit and minFit array to the value of the fit
	for(i = 0; i <= (*ptr_data_size)-1; i++)
	{
		maxFit[i] = darr_extrapolated[i];
		minFit[i] = darr_extrapolated[i];
	}
	
//Uncertainty calculation, varies each parameter by +, - or no uncertainty and verifies if max or min at each altitude
	double tempFit[*ptr_data_size];	
	int iCbaPre, iCbaPost, iAmp, iDec, iOme, iPhi, iSlo;
	double d_tempCbaPre, d_tempCbaPost, d_tempAmp,d_tempTau,d_tempOme,d_tempPhi,d_tempSlo;	
	for (iCbaPre = -1; iCbaPre <= 1; iCbaPre++){
	for (iCbaPost = -1; iCbaPost <= 1; iCbaPost++){
	for (iAmp = -1; iAmp <= 1; iAmp++){
	for (iDec = -1; iDec <= 1; iDec++){
	for (iOme = -1; iOme <= 1; iOme++){
	for (iPhi = -1; iPhi <= 1; iPhi++){
	for (iSlo = -1; iSlo <= 1; iSlo++){
		//Vary the parameters
		d_tempCbaPre = d_CbaPre + result.xerror[1]*iCbaPre*iarr_deviationScale[1];	
		d_tempCbaPost = d_CbaPost + result.xerror[2]*iCbaPost*iarr_deviationScale[2];	
		d_tempAmp = d_amp + result.xerror[4]*iAmp*iarr_deviationScale[4];	
		d_tempTau = d_tau + result.xerror[5]*iDec*iarr_deviationScale[5];	
		d_tempOme = d_ome + result.xerror[6]*iOme*iarr_deviationScale[6];	
		d_tempPhi = d_phi + result.xerror[7]*iPhi*iarr_deviationScale[7];	
		d_tempSlo = d_slo + result.xerror[8]*iSlo*iarr_deviationScale[8];	

		/*Step functions*/
		double darr_stepPre[*ptr_data_size];
		double darr_stepPost[*ptr_data_size];

		for(i = 0; i<*ptr_data_size ;i++){
	
			if((i < d_indBA + index_minFit)|| (d_indBA < 0)){ 		//If the indBA is negative, it means that there is no data post BA
				darr_stepPre[i] = 1;
				darr_stepPost[i] = 0;
			}else{
				darr_stepPre[i] = 0;
				darr_stepPost[i] = 1;
			}
		}

		//Evaluating model
		for (i=0; i<*ptr_data_size; i++) {
			poly=-d_tempAmp*sin(d_tempPhi)+d_tempSlo*darr_timeDif[i];
			s=sin(darr_timeDif[i]*d_tempOme+d_tempPhi);
			e=exp(darr_timeDif[i]*d_tempTau);
			tempFit[i] = (d_tempCbaPre*darr_stepPre[i] + d_tempCbaPost*darr_stepPost[i])*d_x0*(1 - (d_tempAmp*e*s+poly));
		}


		//Verify Max or Min fit		
		for(i = 0; i <= *ptr_data_size-1; i++ )
		{
			if(tempFit[i] > maxFit[i])
			{
				maxFit[i] = tempFit[i];
			}

			if(tempFit[i] < minFit[i])
			{
				minFit[i] = tempFit[i];
			}
		}
		
	}}}}}}}

	//	printf("Uncertainty calculated.\n");

//Return the new max and min uncertainty
	for(i = 0; i < *ptr_data_size; i++ )
	{
		//Return the new uncertainty in terms of 'counts' (NOT extinction)
		ptr_maxUncertainty[i] = ptr_signal[i]-maxFit[i] + d_x0;
		ptr_minUncertainty[i] = ptr_signal[i]-minFit[i] + d_x0;
	}


	//	printf("Flag calculations.\n");

//Flag calculation
	*ptr_chisFlag = 0.0;			//Red. Chi^2 is within spec
	//	printf("Chis: %f, %f \n",*ptr_resChis, *ptr_maxChis);

	if(*ptr_resChis > *ptr_maxChis){
		*ptr_chisFlag = 1.0;		//Flag activation
	}

	*ptr_bumpFlag = 0.0;			//Bump is within spec

//Look for bump in the extrapolation.
//The search algorithm will look for a bump from alt 0 to maxfitz, the definition of index_maxFit DEPENDS on the type of event.
//*ptr_mode: 0 = sunrise, 1 = sunset 
	if(*ptr_mode){	//SS
		for(i = 0; i < *ptr_data_size-index_maxFit+index_BAEnd; i++ )
		{
			float ext = (1.0 - (ptr_residual[*ptr_data_size - i -1])/d_x0);
			if(ext < (*ptr_maxBump)){	//Max bump must be in extinction units!
			//	printf("Bump detected at %d with bump size %f\n", i, ext);
				*ptr_bumpFlag = 1.0;		//Flag activation
			}	
		}
	}else{
		for(i = 0; i < index_minFit + index_BAInit; i++ )
		{
			float ext = (1.0 - (ptr_residual[i])/d_x0);
			if(ext < (*ptr_maxBump)){	//Max bump must be in extinction units!
			//	printf("Bump detected at %d with bump size %f\n", i, ext);
				*ptr_bumpFlag = 1.0;		//Flag activation
			}	
		}


	}

	return 1;	//Success
}



//Auto glue function. IDL gives the option of having the parameters passed in the form argc, argv
//Calling convention, denominated AUTO_GLUE in IDL documentation
int sofieNO8(int argc, void* argv[])
{
	if(argc != 20) return 0;

	double* ptr_signal = (double*) argv[0];
	double* ptr_altitude = (double*) argv[1];
	double* ptr_time = (double*) argv[2];
	char* ptr_mode = (char*) argv[3];
	double* ptr_residual = (double*) argv[4];
	int* ptr_data_size = (int*) argv[5];
	double* ptr_resChis = (double*) argv[6];
	double* ptr_maxChis= (double*) argv[7];
	int* ptr_chisFlag= (int*) argv[8] ;
	double* ptr_Aparameters= (double*) argv[9]; 
	double* ptr_SigmaParameters = (double*) argv[10];
	double* ptr_maxUncertainty= (double*) argv[11];
	double* ptr_minUncertainty= (double*) argv[12];
	int* ptr_bumpFlag= (int*) argv[13];
	double* ptr_maxBump= (double*) argv[14];
	double* ptr_weight=(double*) argv[15];
	double* ptr_minFitz=(double*) argv[16];
	double* ptr_maxFitz=(double*) argv[17];
	double* ptr_startBA=(double*) argv[18];
	double* ptr_endBA=(double*) argv[19];

	return sofieNO8_natural(ptr_signal,
			 	ptr_altitude, 
				ptr_time,
				ptr_mode, 
				ptr_residual, 
				ptr_data_size, 
				ptr_resChis,
				ptr_maxChis,
				ptr_chisFlag,
				ptr_Aparameters,
				ptr_SigmaParameters,
				ptr_maxUncertainty,
				ptr_minUncertainty,
				ptr_bumpFlag,
				ptr_maxBump,
				ptr_weight,
				ptr_minFitz,
				ptr_maxFitz,
				ptr_startBA,
				ptr_endBA);
}




/*************************************************/
//////////////Additional Functions/////////////////
/*************************************************/
//Model to be fitted
int model(int m, int n, double *p, double *deviates, double **derivs, void *ptr_vars)
{
/***************************************************************************************
* This routine models the thermal oscillation in the signal (in terms of counts). The following model definition 
* is used: f(X) = [CpreBA*S(t<tba) + CpostBA*S(t>tba)]{1 - [Amp*e^(X*Tao)*sin(X*Omega + phi) - Amp*sin(phi) + Slope*X]}. Where S(t><tba) is a step function
* at the time of the Balance Adjustment (BA) procedure. The model is then multiplied by a fitted constant pre or post BA.  
*
* Directly called by MPFIT.
*
* Required Input (following convention required by MPFIT.c):
*   m............number of data points 
*   n............number of parameters
*   p............address of parameter array with the form [Amp,Tao,Omega,Phi,Slope]
*   deviates.....address of array where the deviates will be stored
*   derivs.......address of array where derivatives are calculated.
*   ptr_vars.....structure of the form ptr_vars containing the data for the model (independent, dependent variables, error, etc)
*	    .....for more information, consult the documentation of mpfit.c
*
* Output:
*
*   deviates.....deviation from required value is returned, this value will be minimized
*   derivs.......derivative calculation, if not specified, then the algorithm uses a finite difference approximation
*
* Authors: Dr. John McNabb, Hampton University
*	   David Gomez, Hampton University
*************************************************************************************/	


	double x0 = p[0];
	double CbaPre = p[1];
	double CbaPost = p[2];
	double indBA = p[3];
	double amp = p[4];
	double tau = p[5];
	double ome = p[6];
	double phi = p[7];
	double slo = p[8];

	//Extract data from the vars structure
	struct vars_struct *v = (struct vars_struct *) ptr_vars;
	double *x, *y, *ey;
	x = v->x;
	y = v->y;
	ey = v->ey;

	

	/*Step functions*/
	double darr_stepPre[m];
	double darr_stepPost[m];

	int i;
	for(i = 0; i<m ;i++){

		if((i < indBA)||(indBA < 0)){ 		//If the indBA is negative, it means that there is no data post BA
			darr_stepPre[i] = 1;
			darr_stepPost[i] = 0;
		}else{
			darr_stepPre[i] = 0;
			darr_stepPost[i] = 1;
		}
	}

	double poly, s, e;
	/* Compute function deviates */
	for (i=0; i<m; i++) {
		poly=-amp*sin(phi)+slo*x[i];
		s=sin(x[i]*ome+phi);
		e=exp(x[i]*tau);
		deviates[i] = ( y[i] - (CbaPre*darr_stepPre[i] + CbaPost*darr_stepPost[i])*x0*(1 - (amp*e*s+poly)) ) / (ey[i]);
	}

	return 0;
}