//Libraries #include #include #include "mpfit.h" //Downloaded from http://www.physics.wisc.edu/~craigm/idl/cmpfit.html //#include "idl_export.h" //Function Declarations int model(int m, int n, double *p, double *deviates, double **derivs, void *private); int evalModel(int m, int n, double *p, double *eval, void *ptr_vars); //Variable structure (used by mpfit) struct vars_struct{ double *x; //Contains the data used in the fit double *y; //Contains the data used in the fit double *ey; //Contains error of the data used in the fit double *yComplete; //Contains the data outside the fitting region (data below the altitude of minFitz) double *xComplete; //Contains the data outside the fitting region int yComplete_n; //Number of elements of the data outside the fitting region int b_penalize; //If set to 1, the unphysical LVL1 extinction penalty is applied int vars_mode; //Mode of the current event }; //Removal routine int sofieNO12_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 (signal data - fitted oscillation) * is calculated and returned to the address indicated by "ptr_residual" (units will be instrument counts). * The model's decay and frequency are fixed parameters in the model definition. * * Version 11 incorporates improvements that allow an expansion of the fitting region selected according to the * best chi square obtained. In other words, the altitude range for which the best fit was obtained is used. Additionally * a penalty for unphysical residual values is added. If the residual signal (signal data - fitted oscillation) yields results * that are not physically possible (the difference is greater than 0, meaning the residual signal is greater than the * exoatmospheric value expected), then the fit is penalized. By penalization it is meant that the fit deviations are * artificially increased based on how many unphysical datapoints there are and the magnitude of their values. These deviations * are the values that the fitting routine (MPFIT) tries to minimize, therefore MPFIT will modify the parameters to minimize the * amount and magnitude of unphysical residual values. * * UNCERTAINTY: The uncertainty calculations include the error in the fitting parameters, the error due to using a model with * a fixed decay or a model with a restricted decay (mean value +- 3 standard deviations), and the uncertainty due to not expanding * the fit to lower altitudes. This last uncertainty takes into account the lowest possible fit (the lowest possible expanded fit) * that has reasonable reduced chis squared and compares it to the actual fit selected (the one with the best reduced chis square). * * * 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 in the array [X0, CbaPre, CbaPost, indBA, Amp,Tao,Omega,Phi,Slope]. * ..........Look at model for more details. * 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. * ptr_minFitz...........Altitude at which the expanded fit ends. *** 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. *Version 9 - Residual calculation correction, peak threshold lowered and model evaluation as a separate function. *Version 10- Add penalty to negative extinction *Version 10_1- Modification, iteratively look for the best altitude for the fit. *Version 11- Performs both or either enhancements from version 10 and 10_1, according to control variables at the beginning of the file. - Control variables: i_performExpand, i_performPenalty, i_calculateUncertainty *Version 12 - Improve uncertainty calculation to account for extended fits that are also valid. /*******************************/ // printf("Initializing function. \n"); //Control variables //Use these variables to modify the behaviour of the function. //The main advantage of these is to reduce the computation time if one of these is not required. int i_performPenalty = 1; //If 1, endoatmospheric negative extinction will increase the fit deviation and force it to reduce the negative values. int i_performExpand = 1; //If 1, a expanded fit is performed that allows for a better characterization of the fit parameters int i_calculateUncertainty = 1; //If 1, the uncertainty is calculated. //Constants const double d_pi = 3.14159265; const float f_zStep = 0.2; //Altitude step from data point to data point (Km) const float f_peakThreshold = 5.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 = 5; //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. //Variable declaration double i_z0; //Altitude (km) at which the fitting starts, defined with d_minfitz 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 darr_extrapolatedUncertainty[*ptr_data_size]; //Used to save the uncertainty double darr_extrapolatedLowestFit[*ptr_data_size]; //Used to save the lowest fit with acceptable chis2 double d_minFitz = *ptr_minFitz; //Min altitude for the fit (Km) double d_maxFitz = *ptr_maxFitz; //Max altitude for the fit (Km) int b_initBAfound; int b_endBAfound; int counterDatapoints; //Counts for a certain amount of datapoints after the beginning of the BA so that the end of BA is correctly identified //Iteration variables for the minimum fit identification double d_minIterationZ = 100.0; //km This is the minimum altitude that will be verified for the expanded fit algorithm double d_iterationStep = 1.0; //km This is the step between altitudes for the expanded fit algorithm 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); 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; int success_; //If the expanded fit is not going to be performed, set the minimum altitude for the expanded fit d_minIterationZ to the minimum fit limit *ptr_minFitz if(i_performExpand == 0){ d_minIterationZ = d_minFitz; }//End if perform expand //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////////////////// ////////////////////////////////////////////// i_z0 = d_minFitz; //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. This is the reference point. 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). Time differential is the variable used for the fitting. for(i = 0; i < *ptr_data_size; i++ ) { darr_timeDif[i] = (ptr_time[i]) - d_t0; } // printf("Calculated time differential.\n"); //Definition of variables used to save the parameters and results of the best fit. The fit is performed several times, increasing the valid altitude range for the fit and selecting only the best one. double savedChis, savedMinFit, tempChis, savedChisUncertainty; int savedBAinit, savedBAend, savedIndexMinFit, savedIndexMinFitUncertainty; //The BA init and end index with respect to the best fitting range. Also the min fit index. double savedParameters[i_nParameters]; //Array to store the best parameters with the default model double savedParametersUncertainty[i_nParameters]; //Array to store the best parameters with a restricted decay savedChis = 1000000000.0; //An maximum value for the failed events savedChisUncertainty = 1000000000.0; tempChis = 0.0; //Values used to restrict the open decay, taken from the temperature data analysis double temperatureDecay = -0.0327554; double temperatureDecayStddev = 0.00294447; //Saved variables for fits that are also probable, for uncertainty calculations //The code will expand the fit as low as it can and will save the lowest possible fit with an acceptable chis 2. Acceptable refers to a chis2 less than lowestAltitudeChisThreshold. double savedParameters_LowestAltitude[i_nParameters]; double savedLowestValidAltitude = d_maxFitz; int savedIndexMinFitLowestValidAltitude = 0; double lowestAltitudeChisThreshold = 3.0; //The chis^2 limit above which the fit assumes it is an error. This value is used to determine the uncertainty due to expanding the fit. //Saved results for the best fit mp_result savedResult; //This is an MPFIT result structure, where the norm, sigmas, and other parameters are saved. bzero(&savedResult,sizeof(savedResult)); //Initialize the structure savedResult.xerror = perror; //Initialize the structure double parameterValues[] = {0, 1, 1, 0, initialGuess[0], initialGuess[1], initialGuess[2], initialGuess[3], initialGuess[4]}; //Initial guess for the fitting parameters mp_result result; //MPFIT result structure bzero(&result,sizeof(result)); //Initialize the structure /****************************START CYCLES TO EXPAND FIT*******************************/ // This do cycle will gradually expand the fitting region by steps (d_iterationStep). Once the minimum altitude is reached (d_minIterationZ), the cycle will stop // and the best fitting result will be reported. do{ //Get first fitting index for(i = 0; i < *ptr_data_size; i++) { //This next section of code 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 are inverted. This was done to simplify a bit the logic further along. //*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 f_peakThreshold) && !b_initBAfound){ //DIFFERENCE IS OVER THRESHOLD! //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. ****/ /****index_BAInit is for either mode, the first point where a large difference was detected. BAend is the ****/ /****second large difference detected (if there is any second difference detected at all, since not all ****/ /****the portions of the event considered will include the complete BA). **************************************/ //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 int BALength=abs(index_BAEnd - index_BAInit); //Calculate the length of the BA in terms of index if(!b_endBAfound && (*ptr_mode)){BALength=abs(index_BAInit);} //Special case for SS with no end of BA int fitElements_woBA = i_fitElements - BALength -1; //Size of array with no BA 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 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]; if(darr_eyFit[0]==darr_signalFit[i]){ } } }//end else }//end for // printf("BA removal procedure completed.\n"); //The following section of the code is an addition to the vars structure. It includes the signal below our minimal fitting value. In other words //it includes the data from 0 to 140km. This will be used to penalize for any negative extinctions found in this region. //Create variable that contains the entire event after the 140km limit for the fit! int i_eventDataPoints; //If it is a sunset, the event extends from index_Maxfit to the last data point. //If it is a sunrise, the event extends from 0 to index_MinFit //This algorithm is inclusive (will take the last point as part of the array) if(*ptr_mode){ //SS i_eventDataPoints = (*ptr_data_size)- index_maxFit; // printf("Fitting about to start Event data points %i", i_eventDataPoints); }else{ //SR i_eventDataPoints = index_minFit; // printf("Fitting about to start Event data points %i", i_eventDataPoints); } //Set the size of the array holding the entire signal post BA double darr_signalPostZ0[i_eventDataPoints]; double darr_timePostZ0[i_eventDataPoints]; //Assign values to array if(*ptr_mode){ //SS for(i=0; i d_minIterationZ); /********************************************END WHILE****************************************/ //Restore variables to best saved parameters for(i = 0; i0){ //If no complete BA adjustment, leave d_indBA as negative. savedParametersUncertainty[3] += savedIndexMinFitUncertainty; } //Evaluate the model for the uncertainty double temp_parameterValues[] = {savedParametersUncertainty[0], savedParametersUncertainty[1], savedParametersUncertainty[2], savedParametersUncertainty[3], savedParametersUncertainty[4], savedParametersUncertainty[5], savedParametersUncertainty[6], savedParametersUncertainty[7], savedParametersUncertainty[8] }; //vars does not change success_ = evalModel(*ptr_data_size,i_nParameters,temp_parameterValues,darr_extrapolatedUncertainty,&vars); //Extrapolate the fit with acceptable chis2 and lowest fitting range //If no acceptable fit was found, set skip this section. if(savedLowestValidAltitude != d_maxFitz){ //Set the adequate position for the BA if(savedParameters_LowestAltitude[3]>0){ //If no complete BA adjustment, leave d_indBA as negative. savedParameters_LowestAltitude[3] += savedIndexMinFitLowestValidAltitude; } //Evaluate the model for the uncertainty temp_parameterValues[0] = savedParameters_LowestAltitude[0]; temp_parameterValues[1] = savedParameters_LowestAltitude[1]; temp_parameterValues[2] = savedParameters_LowestAltitude[2]; temp_parameterValues[3] = savedParameters_LowestAltitude[3]; temp_parameterValues[4] = savedParameters_LowestAltitude[4]; temp_parameterValues[5] = savedParameters_LowestAltitude[5]; temp_parameterValues[6] = savedParameters_LowestAltitude[6]; temp_parameterValues[7] = savedParameters_LowestAltitude[7]; temp_parameterValues[8] = savedParameters_LowestAltitude[8]; //vars does not change success_ = evalModel(*ptr_data_size,i_nParameters,temp_parameterValues,darr_extrapolatedLowestFit,&vars); }else{ for(i = 0; i<*ptr_data_size ;i++){ darr_extrapolatedLowestFit[i] = darr_extrapolated[i]; //This will cancel this contribution. } }//endif lowest fit with acceptable chis square was found }//endif uncertainty calc /*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)||(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; } } //Residual calculation and model uncertainty calculation for (i=0; i<*ptr_data_size; i++) { // ptr_residual[i] = ptr_signal[i] - darr_extrapolated[i] + (d_CbaPre*darr_stepPre[i] + d_CbaPost*darr_stepPost[i])*d_x0; //The following can be set to flatten the difference between the data before and after the BA. However, it is not used for this version. if(*ptr_mode && d_indBA > 0){ ptr_residual[i] = ptr_signal[i] - darr_extrapolated[i] + (d_CbaPost)*d_x0; //SS }else{ ptr_residual[i] = ptr_signal[i] - darr_extrapolated[i] + (d_CbaPre)*d_x0; //SR } //Calculate uncertainty if one does not fix the decay + the uncertainty due to expanding the fit to lower altitudes if(i_calculateUncertainty){ darr_extrapolatedUncertainty[i] = pow( pow(darr_extrapolatedUncertainty[i] - darr_extrapolated[i],2) + pow(darr_extrapolatedLowestFit[i] - darr_extrapolated[i], 2),0.5); }//End if calculate uncertainty }//End for residual and uncertainty calculation // 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= maxFit[i]) { maxFit[i] = tempFit[i]; //Calculate the uncertainty //ptr_maxUncertainty[i] = ptr_signal[i]-maxFit[i] + (d_tempCbaPre*darr_stepPre[i] + d_tempCbaPost*darr_stepPost[i])*d_x0; //The following can be set to flatten the difference between the data before and after the BA. However, it is not used for this version. if(*ptr_mode && d_indBA > 0){ ptr_minUncertainty[i] = ptr_signal[i]-maxFit[i] + d_tempCbaPost*d_x0; //SS }else{ ptr_minUncertainty[i] = ptr_signal[i]-maxFit[i] + d_tempCbaPre*d_x0; //SR } } if(tempFit[i] <= minFit[i]) { minFit[i] = tempFit[i]; //Calculate the uncertainty //ptr_minUncertainty[i] = ptr_signal[i]-minFit[i] + (d_tempCbaPre*darr_stepPre[i] + d_tempCbaPost*darr_stepPost[i])*d_x0; //The following can be set to flatten the difference between the data before and after the BA. However, it is not used for this version. if(*ptr_mode && d_indBA > 0){ ptr_maxUncertainty[i] = ptr_signal[i]-minFit[i] + d_tempCbaPost*d_x0; //SS }else{ ptr_maxUncertainty[i] = ptr_signal[i]-minFit[i] + d_tempCbaPre*d_x0; //SR } } } }}}}}}} //Add the uncertainty due to fixing/limiting the decay for(i = 0; i < *ptr_data_size; i++ ){ ptr_maxUncertainty[i] = ptr_maxUncertainty[i] + fabs(darr_extrapolatedUncertainty[i]); ptr_minUncertainty[i] = ptr_minUncertainty[i] - fabs(darr_extrapolatedUncertainty[i]); } }//Endif calculate uncertainty // printf("Uncertainty calculated.\n"); // 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 and the BA found. // *ptr_mode: 0 = sunrise, 1 = sunset // recall BAinit = (*ptr_mode && b_endBAfound)?savedBAend:savedBAinit; if(*ptr_mode && b_endBAfound){ //This is a sunset and the endBA was found. for(i = 0; i < (*ptr_data_size)-index_minFit-savedBAend; i++ ) // for(i = 0; i < *ptr_data_size-index_maxFit-savedBAend; i++ ) { int tempIndex = *ptr_data_size - i-1; if(tempIndex < 0){ //Invalid, set flag to indicate issue *ptr_bumpFlag = 1.0; break; } float ext = (1.0 - (ptr_residual[tempIndex])/((d_CbaPre*darr_stepPre[tempIndex] + d_CbaPost*darr_stepPost[tempIndex])*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{ if(*ptr_mode){ //it is a sunset, but no BAend found for(i = 0; i < *ptr_data_size-index_minFit-savedBAinit; i++ ) { int tempIndex = *ptr_data_size - i-1; if(tempIndex < 0){ //Invalid, set flag to indicate issue *ptr_bumpFlag = 1.0; break; } float ext = (1.0 - (ptr_residual[tempIndex])/((d_CbaPre*darr_stepPre[tempIndex] + d_CbaPost*darr_stepPost[tempIndex])*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{ //SR for(i = 0; i < index_minFit + savedBAinit; i++ ) { int tempIndex = i; if(tempIndex < 0){ //Invalid, set flag to indicate issue *ptr_bumpFlag = 1.0; break; } float ext = (1.0 - (ptr_residual[tempIndex])/((d_CbaPre*darr_stepPre[tempIndex] + d_CbaPost*darr_stepPost[tempIndex])*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 } } }//endif ptrmode }//endif ptrmode && b_endBAfound 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 sofieNO12(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 sofieNO12_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(ttba)]{1 - [Amp*e^(X*Tao)*sin(X*Omega + phi) - Amp*sin(phi) + Slope*X]}. Where S(t>x; y = v->y; ey = v->ey; double penalty, penaltyTemp, penaltyOffset, penaltySign, stdDeviation; //penaltySign will set a negative penalty to negative values and a positive penalty to positive ones //Penalty offset defines how much variation is allowed to be considered normal (and variations below this value will not be penalized) //The penalty offset and the standard deviation of the data is set to the value 0.5407 which is the fixed approximated error for the signal. //Further improvements to the algorithm can actually incorporate refinements to these values. penaltyOffset = 0.5407; stdDeviation = 0.5407; penalty = 0.0; int pcount=0; //debug variable, counts the number of data points that would cause negative extinction //Check if penalization should be performed //////////////////////////////////////////// if(v->b_penalize){ //Negative extinction penalty //Extrapolate model of the signal below 140km (yComplete) double *yComplete, *xComplete; int yComplete_n; yComplete = v->yComplete; xComplete = v->xComplete; yComplete_n = v->yComplete_n; double tempModel[yComplete_n]; //Create a temporary 'vars' structure that has the x = xComplete and y = yComplete as primary data for the model evaluation struct vars_struct varsTemp; //Contains the data, x,y and ey (weight) varsTemp.x = xComplete; varsTemp.y = 0; varsTemp.ey = 0; //Extrapolate the fit double p3temp=(v->vars_mode && p[3] > 0)?0:-1; //For SS with a BAend, the CBApost constant must be used below 140km, so we set indBA (p3Temp) to 0. //For all other cases, we want the model to use the CBApre for the entire data below 140km, so we set //it to -1 (no BA flag). double pTemp[] = {p[0], p[1], p[2], p3temp, p[4], p[5], p[6], p[7], p[8]}; tempSuccess_ = evalModel(yComplete_n,n,pTemp,tempModel,&varsTemp); //Verify tempSignal for negative extinction //Compute tempResidual = Signal-Model, if there is a positive difference (residual > 0), it implies that there will be a negative extinction //and then the fit is penalized float penaltyMax=0; int penaltyMaxi=-1; float penaltyRat=1; for (i=0; i penaltyOffset){ //Addition of all values above threshold penalty += penaltyTemp-penaltyOffset; pcount++; //Debug variable } if(penaltyTemp > penaltyMax){penaltyMax=penaltyTemp;penaltyMaxi=i;penaltyRat=yComplete[i]/tempModel[i];} }//EndFOR penalty computation }//EndIF perform penalty? //////////////////////////////////////////// //Compute deviates (mpfit minimizes the deviates value) for (i=0; i= 0){ penaltySign = 1.0; } //Adding the penalty deviates[0] += penalty*penaltySign/(stdDeviation); return 0; } /**************************************************************************************/ int evalModel(int m, int n, double *p, double *eval, void *ptr_vars) { /*************************************************************************************** * This routine evaluates the model. * * Required Input (): * m............number of data points * n............number of parameters * p............address of parameter array with the form [X0, CbaPre, CbaPost, indBA, Amp,Tao,Omega,Phi,Slope] * eval.........address of array where the evaluation results 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: * * eval...array with the resulting evaluation of the model * * 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; int i; x = v->x; y = v->y; ey = v->ey; /*Step functions*/ double darr_stepPre[m]; double darr_stepPost[m]; for(i = 0; i