/********************************************************** * * This program simulates light curves, for 31-parameter (raw - direct from DEBiL program) * * compile: gcc timing.c -o timing -Wall -lm -O4 * * run: timing cmDra2.dbl cmDra2.res 0.4465 0.1855 * **********************************************************/ #include #include #include #include // for strlen(), memmove() and strcpy() #include #define DO_AGC // do auto-gain control (AGC) #define DO_PRINT_OmC_PLOT "plot" // filename extension for plotting results. Used in SM #define DO_PRINT_FOURIER "fourier" #define DO_PRINT_BEST_FIT_SIN "fitSin" #define DO_PRINT_BEST_FIT_PARAB "fitParab" #define NUM_PRINT_BEST_FIT 1000 //#define DEBUG #define OUTLIER_SIGMA_LIMIT 6.0 // for "sigma clipping" #define FOURIER_SPAN_FACTOR 100.0 // should be at least 1.0 #define MAX_FOURIER_FREQ 0.25 // in units of the EB period, to prevent over-fitting #define INTEGRATION_STEP 0.01 // acuracy of light curve simulation //Eccentric anomaly calculation parameter: #define ECC_ANOMALY_MAX_ERROR 1.0e-4 // Convergence threshold (average error is less than a third of this) #define EPSILON 1.0e-8 // to avoid division by zero #define BRACKET_MAX_ERROR 1.0e-6 // must not be so small that there will be the rounding error: (BRACKET_MAX_ERROR+1==1) #define MAX_MAGNITUDE 40.0 // General bisection implementation parameter: // only needed for pathological cases (e.g. overflows or huge roundoff errors) #define BISECTION_SAFETY_MAX_ITERATIONS 1000 #define DEFAULT_LIMB_QUAD_A 0.3542 // 0.282 [Claret 2003 // Claret 1998] #define DEFAULT_LIMB_QUAD_B 0.1939 // 0.311 // A useful constant for magnitude error conversion (2.5 / ln(10)) #define MAGNITUDE_ERROR_FACTOR 1.085736205 #define GOLDEN_RATIO 0.381966011 // = (3 - sqrt(5)) / 2 // Must be at least 2, since otherwise you can always make a fit with chi2-->0, thus having an undefined uncertainty #define MIN_NUM_DATA_IN_DIP 3 #define MIN_NUM_OmC 6 // must be at least 6, for the sinusoidal F-test #define DIP_FIND_STEP_SIZE 0.002 // much smaller than dip "valley", for the initial search #define FRACTION_OF_PLATEAU_TO_USE 0.25 // for auto-gain-control: #define MIN_NUM_DATA_IN_AGC_PROBE 3 #define DO_DIP_IF_AGC_PROBE_FAILED #define DO_ZERO_UNCERTAINTY_IF_MEASUREMENT_FAILED #define MAX_STRING_LEN 256 #define NUM_BOXES_X 1 // when plotting #define NUM_BOXES_Y 4 #define OmC_MULTIPLIER 86400.0 // seconds per day #define BASE_TIME_ROUNDING_FACTOR 100.0 // #define USE_OGLE_LC_FORMAT #define STATISTIC_ERROR_VALUE -100.0 double sqr (double x) { return (x * x) ; } double cube (double x) { return (x * x * x) ; } // Returns the modulo 1 value of x double mod1 (double x) { return (x - floor(x)) ; } // Swaps two given values void swap (double *x1, double *x2) { double tmp = *x1 ; *x1 = *x2 ; *x2 = tmp ; } char* removePath (char *filename) { int i, j ; for (j = 0, i = 0 ; i < strlen(filename) ; i++) if (filename[i] == '/') j = i+1 ; return (&filename[j]) ; } //------------------ // Data structure of a single light curve observation struct LCdata { double primary ; double secondary ; double tertiary ; }; // A comparator for a pair-o-double // Returns: -1 if x.primary < y.primary // 0 if x.primary = y.primary // 1 if x.primary > y.primary int comparLCdata (const void *x, const void *y) { if (((struct LCdata*)x)->primary < ((struct LCdata*)y)->primary) return (-1) ; else return (((struct LCdata*)x)->primary > ((struct LCdata*)y)->primary) ; } // Returns 0 = success ; 1 = failure (malloc) int sortSamples (double *time, double *amp, double *err, int size) { int i ; struct LCdata *arr = (struct LCdata*)malloc (size * sizeof(struct LCdata)) ; if (!arr) return (1) ; for (i = 0 ; i < size ; i++) { arr[i].primary = time[i] ; arr[i].secondary = amp[i] ; arr[i].tertiary = err[i] ; } qsort (arr, size, sizeof(struct LCdata), comparLCdata) ; for (i = 0 ; i < size ; i++) { time[i] = arr[i].primary ; amp[i] = arr[i].secondary ; err[i] = arr[i].tertiary ; } free (arr) ; return (0) ; } //--------------------- // Data structure of a single eclipse O-C struct OmCdata { double time ; // in period units (i.e. number of cycles) double OmC ; // in period units double OmC_errUp ; // size of error (in period units) double OmC_errDown ; // size of error (in period units) double AGCfactor ; int eclipseTag ; // 1 = primary eclipse ; 2 = secondary eclipse }; // A comparator for a pair-o-double // Returns: -1 if x.time < y.time // 0 if x.time = y.time // 1 if x.time > y.time int comparOmCdata (const void *x, const void *y) { if (((struct OmCdata*)x)->time < ((struct OmCdata*)y)->time) return (-1) ; else return (((struct OmCdata*)x)->time > ((struct OmCdata*)y)->time) ; } //-------------------------------------------------------------- // Uses a combined Newton-Raphson + bisection method for finding the root (E) of: // M = E - (e * sin(E)) [Kepler's equation] // Given: M = mean anomaly ; e = eccentricity of orbit (0 <= e < 1) // Returns: E = Eccentric Anomaly double eccentricAnomaly (double M, double e) { int counter = BISECTION_SAFETY_MAX_ITERATIONS ; double Emin = M - 1.0, Emax = M + 1.0, E = M ; double f, dfm, dE ; do { f = M + (e * sin(E)) - E ; // May be optimized by setting cos=sqrt(1-sin^2) dfm = 1.0 - (e * cos(E)) ; // Minus differential of f (always non-negative) dE = dfm * ECC_ANOMALY_MAX_ERROR ; if (f > dE) { Emin = E ; dE = Emax - Emin ; if ((dfm * dE) > f) E += (f / dfm) ; else E = 0.5 * (Emax + Emin) ; } else if (f < (-dE)) { Emax = E ; dE = Emax - Emin ; if ((dfm * dE) > -f) E += (f / dfm) ; else E = 0.5 * (Emax + Emin) ; } else return (E) ; } while ((dE > ECC_ANOMALY_MAX_ERROR) && ((--counter) > 0)) ; return (E) ; // should almost never exits here } //========================== Limb Darkening ==================================== // Quadratic limb darkening coefficients for a solar-like star in I filter (Claret 1998) // Note that a "square root" law give a better fit (especially for I filter), but // but the error in assuming a solar star, dominates by far, and the total CPU requirements // would almost double (doing a linear fit, on the other hand, won't save much CPU). static double limbConstA = 0.0 ; static double limbConstB = 0.0 ; static double limbConstC = 0.0 ; static double limbConstD = 0.0 ; static double limbConstE = 0.0 ; static double limbConstF = 0.0 ; static double limbConstG = 0.0 ; static double limbConstH = 0.0 ; // sets the quadratic limb darkening parameters to arbitrary given values. void setLimbDarkening (double quadA, double quadB) { limbConstA = 1.0 - quadA - quadB ; limbConstB = quadA + quadB + quadB ; limbConstC = -quadB ; limbConstD = M_PI * (1.0 - quadA - quadB - quadB) ; limbConstE = 0.5 * M_PI * quadB ; limbConstF = 2.0 * M_PI * (quadA + quadB + quadB) / 3.0 ; limbConstG = M_PI * (1.0 - quadA - quadB) ; limbConstH = M_PI * (6.0 - quadA - quadA - quadB) / 6.0 ; } // cosAngle2 = the square of the cosine of the angle between the zenith of // the star to the given point on the surface (relative to the star's center) // returns the solar limb darkening coefficient // note that it is normalized so that at the star's zenith (cosAngle2 = 1), it returns 1 double limbDarkening (double cosAngle2) { return (limbConstA + (limbConstB * sqrt(cosAngle2)) + (limbConstC * cosAngle2)) ; } // Analytically integrates the limb darkened disk, from r=0 to r=rBack // rBack = the full radius of the star's disk double integrateWholeDisk (double rBack) { return (limbConstH * rBack * rBack) ; } // Analytically integrates the limb darkened disk, from 0 to the given finish // finish = where to end the integration // rBack = the radius of the back star's disk (the one seen as a crescent) double integrateDiskFromStart (double finish, double rBack) { const double x = sqr(finish / rBack) ; return (rBack * rBack * (((limbConstD + (limbConstE * x)) * x) + (limbConstF * (1.0 - cube(sqrt(1.0 - x)))))) ; } // Analytically integrates the limb darkened disk, from the given start to rBack // start = from where to begin the integration // rBack = the radius of the back star's disk (the one seen as a crescent) double integrateDiskToFinish (double start, double rBack) { const double x = 1.0 - sqr(start / rBack) ; return (rBack * rBack * (((limbConstG - (limbConstE * x)) * x) + (limbConstF * cube(sqrt(x))))) ; } // Makes a simple trapezoid numerical integration with varying step size // (smaller, the closer an asimptote) for finding the brightness of // a partially hidden star. This method works better than Simpson's. // start = where to begin the integration // finish = where to end the integration // rFront = the radius of the front star's disk (the one seen as a disk) // rBack = the radius of the back star's disk (the one seen as a crescent) // D = the distance between the centers of the two stars // dr0 = the infinitesimal integration step scale // Assumes: D + rFront >= finish >= start >= |D - rFront| >= 0 // [num iterations] = 4 / dr0 ; [error] = 0.1 * rFront^2 * dr0^2 double integrateCrescent (double start, double finish, double rFront, double rBack, double D, double dr0) { const double rBack_m2 = 1.0 / (rBack * rBack) ; const double rFront2mD2 = (rFront * rFront) - (D * D) ; const double middle = 0.5 * (start + finish) ; const double D2 = 2.0 * D ; const double limitUp = ((D + rFront < rBack) ? D + rFront : rBack) ; const double limitDown = fabs(D - rFront) ; double dr, r, r2, rStep, sum = 0.0 ; if (D < EPSILON) // do avoid division by zero return (0.0) ; dr0 *= sqrt(rFront) ; // since dr0 is unitless and we want dr to have // units of "distance" (1 = sum semi-major axese) // 1. integrate from middle up dr = dr0 * sqrt(limitUp - middle) ; for (rStep = middle + dr ; rStep < finish ; rStep += dr) { r = rStep - (0.5 * dr) ; r2 = r * r ; sum += (r * acos((rFront2mD2 - r2) / (D2 * r)) * limbDarkening(1.0 - (rBack_m2 * r2)) * dr) ; dr = dr0 * sqrt(limitUp - r) ; } // 2. add sliver at upper edge r = 0.5 * (finish + rStep - dr) ; r2 = r * r ; dr += (finish - rStep) ; sum += (r * acos((rFront2mD2 - r2) / (D2 * r)) * limbDarkening(1.0 - (rBack_m2 * r2)) * dr) ; // 3. integrate from middle down dr = dr0 * sqrt(middle - limitDown) ; for (rStep = middle - dr ; rStep > start ; rStep -= dr) { r = rStep + (0.5 * dr) ; r2 = r * r ; sum += (r * acos((rFront2mD2 - r2) / (D2 * r)) * limbDarkening(1.0 - (rBack_m2 * r2)) * dr) ; dr = dr0 * sqrt(r - limitDown) ; } // 4. add sliver at bottom edge r = 0.5 * (start + rStep + dr) ; r2 = r * r ; dr += (rStep - start) ; sum += (r * acos((rFront2mD2 - r2) / (D2 * r)) * limbDarkening(1.0 - (rBack_m2 * r2)) * dr) ; return (2.0 * sum) ; } // =============== calculate the observed radiation flux from the binary system =============== // rFront = the radius of the star in front // rBack = the radius of the star in back // D = the distance between the centers of the two stars // dr = the infinitesimal integration step // returns the relative brightness of the crescent of the bigger star, // when it is partially eclipsed by the smaller star double integrateBackOverlap (double rFront, double rBack, double D, double dr) { if (rFront < D) { if (rBack > D + rFront) return (integrateDiskFromStart (D - rFront, rBack) + integrateCrescent (D - rFront, D + rFront, rFront, rBack, D, dr) + integrateDiskToFinish (D + rFront, rBack)) ; // E return (integrateDiskFromStart (D - rFront, rBack) + integrateCrescent (D - rFront, rBack, rFront, rBack, D, dr)) ; // B } if (rFront > D + rBack) return (0.0) ; // D if (rBack > D + rFront) return (integrateCrescent (rFront - D, rFront + D, rFront, rBack, D, dr) + integrateDiskToFinish (rFront + D, rBack)) ; // F return (integrateCrescent (rFront - D, rBack, rFront, rBack, D, dr)) ; // C } /****************************************************************** * returns the total expected flux from both stars of the binary * ignores limb darkening, gravity darkening, etc. * time = scan parameter 0 <= time < 1 * p[DIM] = parameter vector * dr = the infinitesimal integration step * doB12max = 1:recalculate the light curve "plateau" brightness ; 0:don't * * NOTE: [definition of 1 unit distance] * the physical distance between the stars at perihelion (minimal) is (1-e) units * thus the physical distance between the stars at aphelion (maximal) is (1+e) units * in other words, 1 unit is the sum of the distances of the two semi-major axes * ******************************************************************/ double flux (double time, double time0, double sin_i, double r1, double r2, const double B1, const double B2, double e, double omega, double dr) { const double meanAnomaly = 2.0 * M_PI * (time - time0) ; double D2, D, E ; E = eccentricAnomaly(meanAnomaly, e) ; D2 = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; D = sqrt (D2) ; if (D >= r1 + r2) // side by side return ((B1 * integrateWholeDisk(r1)) + (B2 * integrateWholeDisk(r2))) ; if (sin(E + omega) > 0.0) // is star 2 in front? return ((B2 * integrateWholeDisk(r2)) + (B1 * integrateBackOverlap(r2, r1, D, dr))) ; // star 1 is in front return ((B1 * integrateWholeDisk(r1)) + (B2 * integrateBackOverlap(r1, r2, D, dr))) ; } //--------------- minimize the theoretical projected binary distance (D^2) --------------- // note: yes, I know that you can optimize the bracketing algorithm by using // golden ratio cuts. I felt that this would make things unnecessarily // complicated with little benefit. double getDipCenter (double initDip, double time0, double sin_i, double e, double omega) { double E, bracketUpX, bracketDownX, bracketMidX, bracketUpY, bracketDownY, bracketMidY ; double testX, testY ; // step 1: measure slope direction bracketDownX = initDip - (0.5 * DIP_FIND_STEP_SIZE) ; E = eccentricAnomaly(2.0 * M_PI * (bracketDownX - time0), e) ; bracketDownY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; bracketUpX = initDip + (0.5 * DIP_FIND_STEP_SIZE) ; E = eccentricAnomaly(2.0 * M_PI * (bracketUpX - time0), e) ; bracketUpY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; // step 2: bracket the minimum if (bracketDownY < bracketUpY) { bracketMidX = bracketUpX ; bracketMidY = bracketUpY ; do { bracketUpX = bracketMidX ; bracketUpY = bracketMidY ; bracketMidX = bracketDownX ; bracketMidY = bracketDownY ; bracketDownX -= DIP_FIND_STEP_SIZE ; E = eccentricAnomaly(2.0 * M_PI * (bracketDownX - time0), e) ; bracketDownY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; } while (bracketDownY < bracketMidY) ; } else // (bracketDownY > bracketUpY) { bracketMidX = bracketDownX ; bracketMidY = bracketDownY ; do { bracketDownX = bracketMidX ; bracketDownY = bracketMidY ; bracketMidX = bracketUpX ; bracketMidY = bracketUpY ; bracketUpX += DIP_FIND_STEP_SIZE ; E = eccentricAnomaly(2.0 * M_PI * (bracketUpX - time0), e) ; bracketUpY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; } while (bracketUpY < bracketMidY) ; } // step 3: find minimum (golden section search) while (((bracketUpY - bracketMidY) > BRACKET_MAX_ERROR) || ((bracketDownY - bracketMidY) > BRACKET_MAX_ERROR)) { testX = bracketMidX + (GOLDEN_RATIO * (bracketDownX - bracketMidX)) ; E = eccentricAnomaly(2.0 * M_PI * (testX - time0), e) ; testY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; if (testY < bracketMidY) { bracketUpX = bracketMidX ; bracketUpY = bracketMidY ; bracketMidX = testX ; bracketMidY = testY ; } else { bracketDownX = testX ; bracketDownY = testY ; } testX = bracketMidX + (GOLDEN_RATIO * (bracketUpX - bracketMidX)) ; E = eccentricAnomaly(2.0 * M_PI * (testX - time0), e) ; testY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; if (testY < bracketMidY) { bracketDownX = bracketMidX ; bracketDownY = bracketMidY ; bracketMidX = testX ; bracketMidY = testY ; } else { bracketUpX = testX ; bracketUpY = testY ; } } return (bracketMidX) ; } double getDipHalfWidth (double dipMinimum, double stepSize, double time0, double sin_i, double e, double omega, double r1r2) { double bracketDownX, bracketUpX, bracketDownY, bracketUpY, testX, testY, E ; // step 1: bracket edge of dip r1r2 = sqr(r1r2) ; bracketUpX = dipMinimum ; //------------ do // not that we need to protect against r1r2 ~= 1, for which the bracket might miss { E = eccentricAnomaly(2.0 * M_PI * (bracketUpX - time0), e) ; bracketUpY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; do { bracketDownX = bracketUpX ; bracketDownY = bracketUpY ; bracketUpX += stepSize ; E = eccentricAnomaly(2.0 * M_PI * (bracketUpX - time0), e) ; bracketUpY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; } while ((bracketUpY < r1r2) && (bracketDownY < bracketUpY)) ; if (bracketUpY < r1r2) { if (stepSize < BRACKET_MAX_ERROR) return (fabs(dipMinimum - bracketDownX)) ; bracketUpX -= stepSize ; // two steps back for the worst case scenerio bracketUpX -= stepSize ; stepSize *= 0.2 ; } } while (bracketUpY < r1r2) ; // step 2: find edge while (fabs(bracketDownX - bracketUpX) > BRACKET_MAX_ERROR) { testX = 0.5 * (bracketDownX + bracketUpX) ; E = eccentricAnomaly(2.0 * M_PI * (testX - time0), e) ; testY = sqr(1.0 - (e * cos(E))) - sqr((((cos(E) - e) * sin(omega)) + (sqrt(1.0 - (e * e)) * sin(E) * cos(omega))) * sin_i) ; if (testY > r1r2) bracketUpX = testX ; else bracketDownX = testX ; } return (fabs(dipMinimum - bracketDownX)) ; } // ----------------- minimize the chi2 function (fit to data) -------------------- // scans LC point from indexStart to (indexFinish-1), inclusive double getChi2AGC (int indexStart, int indexFinish, double *time, double *amp, double *err, double time0, double sin_i, double r1, double r2, double B1, double B2, double e, double omega, double AGCfactor) { int i ; double sum = 0.0 ; for (i = indexStart ; i < indexFinish ; i++) sum += sqr(((AGCfactor * amp[i]) - flux (time[i], time0, sin_i, r1, r2, B1, B2, e, omega, INTEGRATION_STEP)) / err[i]) ; return (sum) ; } // assumes that the *time array is sorted // returns the index that is equal or just above the goatTime int getIndexUp (double goalTime, double *time, int size) { int testIndex, bracketDownIndex = 0, bracketUpIndex = size - 1 ; while ((bracketUpIndex - bracketDownIndex) > 1) { testIndex = (bracketUpIndex + bracketDownIndex) / 2 ; if (time[testIndex] < goalTime) bracketDownIndex = testIndex ; else bracketUpIndex = testIndex ; } return (bracketUpIndex) ; } // A comparator for a pair-o-double // Returns: -1 if x < y // 0 if x = y // 1 if x > y int comparDouble (const void *x, const void *y) { if (*((double*)x) < *((double*)y)) return (-1) ; else return (*((double*)x) > *((double*)y)) ; } // return an exact value for dip center // returns MAXFLOAT upon an error double getLocalDipShift (double tDip, double tWidthDown, double tWidthUp, double *time, double *amp, double *err, double *arrPlateau, int size, double time0, double sin_i, double r1, double r2, double B1, double B2, double e, double omega, double lenPlateauPrev, double lenPlateauPost, double *pUncertaintyUp, double *pUncertaintyDown, double *pAGCfactor) { int indexStart, indexFinish, indexStartPlateau, indexFinishPlateau ; double bracketDownX, bracketDownY, bracketUpX, bracketUpY, bracketMidX, bracketMidY ; double testX, testY, markDownX, markUpX, chi2_1sigma ; #ifdef DO_AGC int i, j, numIndexPlateau ; double ampPlateau ; #endif // step 1: get dip index range indexStart = getIndexUp (tDip - tWidthDown, time, size) ; indexFinish = getIndexUp (tDip + tWidthUp, time, size) ; if (indexFinish - indexStart < MIN_NUM_DATA_IN_DIP) return (MAXFLOAT) ; // step 2: get plateau index range indexStartPlateau = getIndexUp (tDip - tWidthDown - lenPlateauPrev, time, size) ; indexFinishPlateau = getIndexUp (tDip + tWidthUp + lenPlateauPost, time, size) ; #ifdef DO_AGC numIndexPlateau = indexFinishPlateau - indexFinish + indexStart - indexStartPlateau - 1 ; if (numIndexPlateau < MIN_NUM_DATA_IN_AGC_PROBE) { #ifdef DO_DIP_IF_AGC_PROBE_FAILED #ifdef DEBUG printf ("Warning: Not enough data points in the eclipse plateau probe to locate it (%d), using AGC=1.0\n", numIndexPlateau) ; #endif *pAGCfactor = 1.0 ; #else #ifdef DEBUG printf ("ERROR: Not enough data points in the eclipse plateau probe to locate it (%d)\n", numIndexPlateau) ; #endif return (MAXFLOAT) ; #endif } else // the AGC plateau probe is OK { // step 3: find AGC from plateau median j = 0 ; for (i = indexStartPlateau ; i < indexStart ; i++, j++) arrPlateau[j] = amp[i] ; for (i = indexFinish + 1 ; i < indexFinishPlateau ; i++, j++) arrPlateau[j] = amp[i] ; qsort (arrPlateau, numIndexPlateau, sizeof(double), comparDouble) ; if (numIndexPlateau & 1) // is odd ampPlateau = arrPlateau[numIndexPlateau / 2] ; else // is even ampPlateau = 0.5 * (arrPlateau[numIndexPlateau / 2] + arrPlateau[(numIndexPlateau / 2) - 1]) ; *pAGCfactor = ((B1 * integrateWholeDisk(r1)) + (B2 * integrateWholeDisk(r2))) / ampPlateau ; #ifdef DEBUG printf ("numIndexPlateau = %d ampPlateau = %e AGC = %f (%d)\n", numIndexPlateau, ampPlateau, *pAGCfactor, numIndexPlateau) ; #endif } #else *pAGCfactor = 1.0 ; #endif // step 4: get dip bracket: bracketDownX = time0 - tWidthDown ; bracketDownY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, bracketDownX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; bracketMidX = time0 ; bracketMidY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, bracketMidX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; bracketUpX = time0 + tWidthUp ; bracketUpY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, bracketUpX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; if ((bracketDownY < bracketMidY) || (bracketUpY < bracketMidY)) { #ifdef DEBUG printf ("ERROR: The bracket fit at the edge is better than the middle (time=%f) [%e %e %e]\n", tDip, bracketDownY, bracketMidY, bracketUpY) ; #endif return (MAXFLOAT) ; } // step 5: find minimum (golden section search) while (((bracketUpY - bracketMidY) > BRACKET_MAX_ERROR) || ((bracketDownY - bracketMidY) > BRACKET_MAX_ERROR)) { testX = bracketMidX + (GOLDEN_RATIO * (bracketDownX - bracketMidX)) ; testY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, testX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; if (testY < bracketMidY) { bracketUpX = bracketMidX ; bracketUpY = bracketMidY ; bracketMidX = testX ; bracketMidY = testY ; } else { bracketDownX = testX ; bracketDownY = testY ; } testX = bracketMidX + (GOLDEN_RATIO * (bracketUpX - bracketMidX)) ; testY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, testX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; if (testY < bracketMidY) { bracketDownX = bracketMidX ; bracketDownY = bracketMidY ; bracketMidX = testX ; bracketMidY = testY ; } else { bracketUpX = testX ; bracketUpY = testY ; } } // this normalizes the minimum chi2 to N. Formally, it should be: bracketMidY + 1.0 chi2_1sigma = bracketMidY + (bracketMidY / (indexFinishPlateau - indexStartPlateau)) ; // step 6a: find width at min+const (down side) bracketUpX = bracketMidX ; bracketDownX = bracketMidX - tWidthDown ; testY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, bracketDownX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; if (testY < chi2_1sigma) { #ifdef DO_ZERO_UNCERTAINTY_IF_MEASUREMENT_FAILED markDownX = bracketDownX ; #else #ifdef DEBUG printf ("ERROR: Finding the uncertainty (down)\n") ; #endif return (MAXFLOAT) ; #endif } else { while ((bracketUpX - bracketDownX) > BRACKET_MAX_ERROR) { testX = 0.5 * (bracketDownX + bracketUpX) ; testY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, testX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; if (testY < chi2_1sigma) bracketUpX = testX ; else bracketDownX = testX ; } markDownX = 0.5 * (bracketDownX + bracketUpX) ; } // step 6b: find width at min+const (up side) bracketUpX = bracketMidX + tWidthUp ; bracketDownX = bracketMidX ; testY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, bracketUpX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; if (testY < chi2_1sigma) { #ifdef DO_ZERO_UNCERTAINTY_IF_MEASUREMENT_FAILED markUpX = bracketUpX ; #else #ifdef DEBUG printf ("ERROR: Finding the uncertainty (up)\n") ; #endif return (MAXFLOAT) ; #endif } else { while ((bracketUpX - bracketDownX) > BRACKET_MAX_ERROR) { testX = 0.5 * (bracketDownX + bracketUpX) ; testY = getChi2AGC (indexStartPlateau, indexFinishPlateau, time, amp, err, testX, sin_i, r1, r2, B1, B2, e, omega, *pAGCfactor) ; if (testY < chi2_1sigma) bracketDownX = testX ; else bracketUpX = testX ; } markUpX = 0.5 * (bracketDownX + bracketUpX) ; } //--------------- *pUncertaintyUp = markUpX - bracketMidX ; *pUncertaintyDown = bracketMidX - markDownX ; return (bracketMidX - time0) ; } //================================================== #ifdef DO_PRINT_OmC_PLOT void writeSM (FILE *foutSM, char *filenameSM, char *filenamePlotBase, char *filenamePlot1, char *filenamePlot2, double period, long baseTime, int hasPrimaryEclipses, int hasSecondaryEclipses, char *filenameBestFitParab, char *filenameBestFitSin) { static int isFirst = 1, i = 0, j = NUM_BOXES_Y ; if (!hasPrimaryEclipses && !hasSecondaryEclipses) return ; if (isFirst) { isFirst = 0 ; fprintf (foutSM, "device postportfile %s.ps\n", filenameSM) ; // fprintf (foutSM, "expand 0.4\n") ; fprintf (foutSM, "ticksize 0 0 0 0\n") ; fprintf (foutSM, "ctype black\n\n") ; fprintf (foutSM, "ltype 0\n") ; } i++ ; if (i > NUM_BOXES_X) { i = 1 ; j-- ; } if (j <= 0) { j = NUM_BOXES_Y ; fprintf (foutSM, "page\n\n") ; } fprintf (foutSM, "window %d %d %d %d\n", NUM_BOXES_X, NUM_BOXES_Y, i, j) ; fprintf (foutSM, "set x = 0\n") ; fprintf (foutSM, "set y = 0\n") ; fprintf (foutSM, "set x1 = 0\n") ; fprintf (foutSM, "set y1 = 0\n") ; fprintf (foutSM, "set y1up_err = 0\n") ; fprintf (foutSM, "set y1down_err = 0\n") ; fprintf (foutSM, "set x2 = 0\n") ; fprintf (foutSM, "set y2 = 0\n") ; fprintf (foutSM, "set y2up_err = 0\n") ; fprintf (foutSM, "set y2down_err = 0\n") ; if (hasPrimaryEclipses) { fprintf (foutSM, "data %s\n", filenamePlot1) ; fprintf (foutSM, "read {x1 1 y1 2 y1up_err 3 y1down_err 4}\n") ; fprintf (foutSM, "set base_x1 = x1 - %ld\n", baseTime) ; } if (hasSecondaryEclipses) { fprintf (foutSM, "data %s\n", filenamePlot2) ; fprintf (foutSM, "read {x2 1 y2 2 y2up_err 3 y2down_err 4}\n") ; fprintf (foutSM, "set base_x2 = x2 - %ld\n", baseTime) ; } if (hasPrimaryEclipses && hasSecondaryEclipses) { fprintf (foutSM, "set x = base_x1 concat base_x2\n") ; // fprintf (foutSM, "set y = (y1+y1up_err) concat (y1-y1down_err) concat (y2+y2up_err) concat (y2-y2down_err)\n") ; fprintf (foutSM, "set y = y1 concat y2\n") ; fprintf (foutSM, "limits x y\n") ; } else if (hasPrimaryEclipses) fprintf (foutSM, "limits base_x1 y1\n") ; else // if (hasSecondaryEclipses) fprintf (foutSM, "limits base_x2 y2\n") ; if (hasPrimaryEclipses) { fprintf (foutSM, "ptype 3 3\n") ; fprintf (foutSM, "points base_x1 y1\n") ; fprintf (foutSM, "errorbar base_x1 y1 y1up_err 2\n") ; fprintf (foutSM, "errorbar base_x1 y1 y1down_err 4\n") ; } if (hasSecondaryEclipses) { fprintf (foutSM, "ptype 4 2\n") ; fprintf (foutSM, "points base_x2 y2\n") ; fprintf (foutSM, "errorbar base_x2 y2 y2up_err 2\n") ; fprintf (foutSM, "errorbar base_x2 y2 y2down_err 4\n") ; } if (*filenameBestFitParab) { fprintf (foutSM, "ltype 1\n") ; fprintf (foutSM, "ctype blue\n") ; fprintf (foutSM, "data %s\n", filenameBestFitParab) ; fprintf (foutSM, "read {x 1 y 2}\n") ; fprintf (foutSM, "set base_x = x - %ld\n", baseTime) ; fprintf (foutSM, "connect base_x y\n") ; fprintf (foutSM, "ltype 0\n") ; fprintf (foutSM, "ctype black\n") ; } if (*filenameBestFitSin) { fprintf (foutSM, "ctype red\n") ; fprintf (foutSM, "data %s\n", filenameBestFitSin) ; fprintf (foutSM, "read {x 1 y 2}\n") ; fprintf (foutSM, "set base_x = x - %ld\n", baseTime) ; fprintf (foutSM, "connect base_x y\n") ; fprintf (foutSM, "ctype black\n") ; } if (baseTime > 0) fprintf (foutSM, "xlabel %s Time - %ld (P=%.3f days)\n", filenamePlotBase, baseTime, period) ; else if (baseTime < 0) fprintf (foutSM, "xlabel %s Time + %ld (P=%.3f days)\n", filenamePlotBase, -baseTime, period) ; else // (baseTime == 0) fprintf (foutSM, "xlabel %s Time (P=%.3f days)\n", filenamePlotBase, period) ; // fprintf (foutSM, "ylabel O-C\n") ; fprintf (foutSM, "box\n\n") ; fflush(foutSM) ; } #endif //------------------- sigma-clipping ------------------------ double getSigma (int index, struct OmCdata *OmCarr) { if (OmCarr[index].OmC > 0.0) return (OmCarr[index].OmC / OmCarr[index].OmC_errUp) ; return (OmCarr[index].OmC / OmCarr[index].OmC_errDown) ; } // Removes the outliers (from a weighted first-order regression) from the input data array, // and normalizes the O-C curve (de-trends), so that its average and slope become zero. // *pAlarmScore = the alarm score [Mazeh, Tamuz, and North (2005)] arournd the best-fit linear regression // Returns the final reduced chi-squared (chi2nu) around the linear regression // Upon error returns -1.0 (too few data points, or infinite slope) double sigmaClipping (struct OmCdata *OmCarr, int *pOmCsize, double *pslope, double *pAlarmScore) { int i, doClip, numClipped = 0, OmCsize = *pOmCsize ; double N, A, B, C, D, E, tmp, slope, offset, chi2nu, t, x, invSigma2, tAvr ; double sumRun, startVal, val, sumAlarm, S2, sigma2, diff ; *pAlarmScore = STATISTIC_ERROR_VALUE ; // should never be outputted *pslope = STATISTIC_ERROR_VALUE ; do { if (OmCsize < MIN_NUM_OmC) { printf ("Orig = %d Num clipped = %d\n", *pOmCsize, numClipped) ; return (-1.0) ; } doClip = 0 ; N = 0.0 ; A = 0.0 ; B = 0.0 ; C = 0.0 ; D = 0.0 ; E = 0.0 ; tAvr = 0.0 ; for (i = 0 ; i < OmCsize ; i++) // for better numerical accuracy tAvr += OmCarr[i].time ; tAvr /= OmCsize ; for (i = 0 ; i < OmCsize ; i++) { t = OmCarr[i].time - tAvr ; x = OmCarr[i].OmC ; invSigma2 = 1.0 / sqr(OmCarr[i].OmC_errUp + OmCarr[i].OmC_errDown) ; // an approximation to a (symmertric) Gaussian N += invSigma2 ; A += (t * invSigma2) ; B += (t * t * invSigma2) ; C += (t * x * invSigma2) ; D += (x * invSigma2) ; E += (x * x * invSigma2) ; } tmp = (N * B) - (A * A) ; if ((tmp < EPSILON) && (tmp > -EPSILON)) return (-1.0) ; slope = ((N * C) - (D * A)) / tmp ; offset = ((D * B) - (C * A)) / tmp ; offset -= (slope * tAvr) ; chi2nu = E - (D * D / N) - (sqr((N * C) - (D * A)) / (tmp * N)) ; chi2nu /= (OmCsize-2) ; #ifdef DEBUG printf ("Num= %d (N=%e) slope= %lf offset= %lf chi2nu= %lf\n", OmCsize, N, slope, offset, chi2nu) ; #endif //---------- for (i = 0 ; i < OmCsize ; i++) { diff = (OmCarr[i].time * slope) + offset - OmCarr[i].OmC ; if (diff > 0.0) sigma2 = sqr(diff / OmCarr[i].OmC_errUp) ; else sigma2 = sqr(diff / OmCarr[i].OmC_errDown) ; if (sigma2 > (OUTLIER_SIGMA_LIMIT * OUTLIER_SIGMA_LIMIT * chi2nu)) { memmove (&OmCarr[i], &OmCarr[i+1], (OmCsize - i - 1) * sizeof(struct OmCdata)) ; OmCsize-- ; i-- ; doClip = 1 ; numClipped++ ; } } } while (doClip) ; printf ("Orig = %d Num clipped = %d\n", *pOmCsize, numClipped) ; *pOmCsize = OmCsize ; *pslope = slope ; // --------- Normalize the O-C values ---------------------- for (i = 0 ; i < OmCsize ; i++) OmCarr[i].OmC -= ((OmCarr[i].time * slope) + offset) ; qsort (OmCarr, OmCsize, sizeof(struct OmCdata), comparOmCdata) ; //------------------------------ calculate the alarm score --------------------- sumAlarm = 0.0 ; i = 0 ; val = getSigma(i, OmCarr) ; while (i < OmCsize) { startVal = val ; sumRun = startVal ; while ((i < OmCsize) && (startVal * (val = getSigma(i, OmCarr)) > 0.0)) // same sign { sumRun += val ; i++ ; } sumAlarm += (sumRun * sumRun) ; } S2 = 0.0 ; for (i = 0 ; i < OmCsize ; i++) S2 += sqr(getSigma(i, OmCarr)) ; *pAlarmScore = (sumAlarm / S2) - 1.0 - (4.0 / M_PI) ; return (chi2nu) ; } //----------- F-test statistics ----------------------------- // Makes a weighted second-order regression (parabola-fit). // Returns the ratio between the linear chi2 and the parabola chi2 (F-test) // *pPdot is the period change over time (in units of seconds per second) // returns zero upon error (when the points are aligned vertically) double calcParabola_Ftest (char *filenameNoPath, double period, struct OmCdata *OmCarr, int OmCsize, double *pPdot, char *filenameBestFitParab) { int i ; double A0 = 0.0, A1 = 0.0, A2 = 0.0, A3 = 0.0, A4 = 0.0 ; double B0 = 0.0, B1 = 0.0, B2 = 0.0, C0 = 0.0 ; double denom, a, b, c, t, t2, x, invSigma2, chi2parabola, tAvr = 0.0 ; #ifdef DO_PRINT_BEST_FIT_PARAB double time ; FILE *foutParab ; #endif *pPdot = STATISTIC_ERROR_VALUE ; // should never be outputted if (OmCsize <= 3) // for F-test return (0.0) ; for (i = 0 ; i < OmCsize ; i++) // for better numerical accuracy tAvr += OmCarr[i].time ; tAvr /= OmCsize ; for (i = 0 ; i < OmCsize ; i++) { t = OmCarr[i].time - tAvr ; t2 = t * t ; x = OmCarr[i].OmC ; invSigma2 = 1.0 / sqr(OmCarr[i].OmC_errUp + OmCarr[i].OmC_errDown) ; // an approximation to a (symmertric) Gaussian A0 += invSigma2 ; A1 += (t * invSigma2) ; A2 += (t2 * invSigma2) ; A3 += (t * t2 * invSigma2) ; A4 += (t2 * t2 * invSigma2) ; B0 += (x * invSigma2) ; B1 += (x * t * invSigma2) ; B2 += (x * t2 * invSigma2) ; C0 += (x * x * invSigma2) ; } denom = (A2 * ((A0*A4) + (2.0*A1*A3) - (A2*A2))) - (A0*A3*A3) - (A1*A1*A4) ; // Denominator if ((denom < EPSILON) && (denom > -EPSILON)) return (0.0) ; // error a = ((A0 * ((A2*B2) - (A3*B1))) + (A1 * ((A3*B0) - (A1*B2))) + (A2 * ((A1*B1) - (A2*B0)))) / denom ; b = ((A0 * ((A4*B1) - (A3*B2))) + (A1 * ((A2*B2) - (A4*B0))) + (A2 * ((A3*B0) - (A2*B1)))) / denom ; c = ((A1 * ((A3*B2) - (A4*B1))) + (A2 * ((A4*B0) - (A2*B2))) + (A3 * ((A2*B1) - (A3*B0)))) / denom ; *pPdot = 2.0 * a ; // Pdot = 2 * P0 * a, but we've normalized: P0 = 1 chi2parabola = (A4*a*a) + (2.0*a*b*A3) + (A2 * ((b*b) + (2.0*a*c))) + (2.0*b*c*A1) + (A0*c*c) - (2.0 * ((B2*a) + (B1*b) + (B0*c))) + C0 ; if ((chi2parabola < EPSILON) && (chi2parabola > -EPSILON)) return (0.0) ; // error #ifdef DEBUG printf ("chi2nu parab= %lf chi2nu linear= %lf F-test= %lf Pdot= %e\n", chi2parabola / (OmCsize-3), C0 / (OmCsize-2), (C0 / (OmCsize-2)) / (chi2parabola / (OmCsize-3)), 2.0 * a) ; #endif //------ print best fit -------------- #ifdef DO_PRINT_BEST_FIT_PARAB sprintf (filenameBestFitParab, "%s_P%.3f.%s", filenameNoPath, period, DO_PRINT_BEST_FIT_PARAB) ; foutParab = fopen (filenameBestFitParab, "wt") ; if (!foutParab) { printf ("ERROR: Couln't open '%s' for writing best parabolic fit\n", filenameBestFitParab) ; *filenameBestFitParab = 0 ; } else { for (i = 0 ; i < NUM_PRINT_BEST_FIT ; i++) { time = OmCarr[0].time + ((OmCarr[OmCsize-1].time - OmCarr[0].time) * (double)i / (NUM_PRINT_BEST_FIT - 1)) ; fprintf (foutParab, "%lf %lf\n", time * period, OmC_MULTIPLIER * period * ((a * sqr(time - tAvr)) + (b * (time - tAvr)) + c)) ; } fclose (foutParab) ; } #else *filenameBestFitParab = 0 ; #endif return ((C0 / (OmCsize-2)) / (chi2parabola / (OmCsize-3))) ; } // Minimizes SUM_i ((factorSin * SIN(w*t_i)) + (factorCos * COS(w*t_i)) - y_i)^2 // If enabled, will write a file for plotting the power spectrum // Note:assumes that the time stamps of OmCarr[] are sorted // returns the best-sinusoidal F-test result (chi2_linear / chi2_sinusoidal) // returns 0, upon error double calcSinusoidal_Ftest (char *filenameNoPath, double period, struct OmCdata *OmCarr, int OmCsize, double *pBestFreq, double *pBestFreqAmp, char *filenameBestFitSin) { double f, w, invSigma2, minFreq, maxFreq, T ; double bestFreq, bestFreqChi2, bestFreqPwr, chi2sin ; double tmp, A, B1, B2, C, D, E, sin_wt, cos_wt, y, factorSin, factorCos ; int i ; FILE *foutFourier = 0 ; #ifdef DO_PRINT_FOURIER char filenameFourier[MAX_STRING_LEN] ; #endif #ifdef DO_PRINT_BEST_FIT_SIN double time, bestFactorSin = 0.0, bestFactorCos = 0.0 ; #endif *pBestFreq = STATISTIC_ERROR_VALUE ; // should never be outputted *pBestFreqAmp = STATISTIC_ERROR_VALUE ; if (OmCsize <= 5) // for F-test return (0.0) ; T = OmCarr[OmCsize-1].time - OmCarr[0].time ; if (T <= 0.0) { printf ("ERROR: Unsorted timestamps\n") ; return (0.0) ; } #ifdef DO_PRINT_FOURIER sprintf (filenameFourier, "%s_P%.3f.%s", filenameNoPath, period, DO_PRINT_FOURIER) ; foutFourier = fopen (filenameFourier, "wt") ; if (!(foutFourier = fopen (filenameFourier, "wt"))) { printf ("ERROR: Couldn't open the Fourier output file\n") ; return (0.0) ; } #endif E = 0.0 ; for (i = 0 ; i < OmCsize ; i++) { y = OmCarr[i].OmC ; invSigma2 = 1.0 / sqr(OmCarr[i].OmC_errUp + OmCarr[i].OmC_errDown) ; // an approximation to a (symmertric) Gaussian E += (y * y * invSigma2) ; } bestFreq = 0.0 ; bestFreqChi2 = E ; bestFreqPwr = 0.0 ; minFreq = 1.0 / (FOURIER_SPAN_FACTOR * T) ; // If it were equal sampling, Nyquist would limit maxFreq = 0.5 * OmCsize / T, beyond which you'd have aliasing. maxFreq = MAX_FOURIER_FREQ ; // to prevent over fitting for (f = minFreq ; f < maxFreq ; f += minFreq) { w = 2.0 * M_PI * f ; A = 0.0 ; B1 = 0.0 ; B2 = 0.0 ; C = 0.0 ; D = 0.0 ; for (i = 0 ; i < OmCsize ; i++) { sin_wt = sin(w * OmCarr[i].time) ; cos_wt = cos(w * OmCarr[i].time) ; y = OmCarr[i].OmC ; invSigma2 = 1.0 / sqr(OmCarr[i].OmC_errUp + OmCarr[i].OmC_errDown) ; // an approximation to a (symmertric) Gaussian A += (sin_wt * cos_wt * invSigma2) ; B1 += (sin_wt * sin_wt * invSigma2) ; B2 += (cos_wt * cos_wt * invSigma2) ; C += (y * sin_wt * invSigma2) ; D += (y * cos_wt * invSigma2) ; } tmp = (B1 * B2) - (A * A) ; if ((tmp < EPSILON) && (tmp > -EPSILON)) return (0.0) ; factorSin = ((B2 * C) - (D * A)) / tmp ; factorCos = ((D * B1) - (C * A)) / tmp ; chi2sin = E - (D * D / B2) - (sqr((B2 * C) - (D * A)) / (tmp * B2)) ; if (foutFourier) // period is used the make the time units into days fprintf (foutFourier, "%lf %lf %lf\n", period / f, f / period, (E / (OmCsize-2)) / (chi2sin / (OmCsize-5))) ; if (chi2sin < bestFreqChi2) { bestFreqChi2 = chi2sin ; bestFreq = f ; bestFreqPwr = (factorSin * factorSin) + (factorCos * factorCos) ; #ifdef DO_PRINT_BEST_FIT_SIN bestFactorSin = factorSin ; bestFactorCos = factorCos ; #endif } } *pBestFreq = bestFreq / period ; *pBestFreqAmp = sqrt(bestFreqPwr) ; if (foutFourier) fclose (foutFourier) ; //------ print best fit -------------- #ifdef DO_PRINT_BEST_FIT_SIN sprintf (filenameBestFitSin, "%s_P%.3f.%s", filenameNoPath, period, DO_PRINT_BEST_FIT_SIN) ; foutFourier = fopen (filenameBestFitSin, "wt") ; if (!foutFourier) { printf ("ERROR: Couln't open '%s' for writing best sinusoidal fit\n", filenameBestFitSin) ; *filenameBestFitSin = 0 ; } else { for (i = 0 ; i < NUM_PRINT_BEST_FIT ; i++) { time = OmCarr[0].time + (T * (double)i / (NUM_PRINT_BEST_FIT - 1)) ; fprintf (foutFourier, "%lf %lf\n", time * period, OmC_MULTIPLIER * period * ((bestFactorSin * sin(2.0 * M_PI * bestFreq * time)) + (bestFactorCos * cos(2.0 * M_PI * bestFreq * time)))) ; } fclose (foutFourier) ; } #else *filenameBestFitSin = 0 ; #endif return ((E / (OmCsize-2)) / (bestFreqChi2 / (OmCsize-5))) ; // 5 fourier parameters: slope, offset, freq, factorSin, factorCos } int main (int argc, char **argv) { FILE *fin, *finLC, *fout ; double quadA = DEFAULT_LIMB_QUAD_A, quadB = DEFAULT_LIMB_QUAD_B ; int size, i, OmCsize, hasPrimaryEclipses, hasSecondaryEclipses ; double period, e, r1, r2, B1, B2, sin_i, time0, omega ; char filename[MAX_STRING_LEN], *filenameNoPath ; char filenameBestFitParab[MAX_STRING_LEN], filenameBestFitSin[MAX_STRING_LEN] ; double tmpTime, tmpMag, tmpErr ; double *time, *amp, *err, *tmpArr ; char formatStr[MAX_STRING_LEN] ; double tDip1, tDip2, tDip, tDip1Wup, tDip1Wdown, tDip2Wup, tDip2Wdown ; double lenPlateau12, lenPlateau21, deltaLocalDip ; double uncertaintyUp, uncertaintyDown, AGCfactor, alarmScore, slope ; double parabFtest, Pdot, chi2nu_linear, sinusoidalFtest, bestFreq, bestFreqAmp ; struct OmCdata *OmCarr ; #ifdef DO_PRINT_OmC_PLOT long baseTime ; char filenamePlot1[MAX_STRING_LEN], filenamePlot2[MAX_STRING_LEN], filenameSM[MAX_STRING_LEN] ; FILE *foutPlot1, *foutPlot2, *foutSM ; #endif if ((argc != 3) && (argc != 5)) { printf ("Usage: %s [ ]\n", argv[0]) ; return (1) ; } if (!(fin = fopen (argv[1], "rt"))) { printf ("ERROR: Couldn't open the DEBiL database\n") ; return (2) ; } if (!(fout = fopen (argv[2], "wt"))) { printf ("ERROR: Couldn't open the output results file\n") ; return (3) ; } if (argc == 5) { quadA = atof(argv[3]) ; quadB = atof(argv[4]) ; } printf ("Limb darkening coefficients: A=%f ; B=%f\n", quadA, quadB) ; setLimbDarkening (quadA, quadB) ; #ifdef USE_OGLE_LC_FORMAT strcpy (formatStr, "%lf %*f %*f %lf %lf %*d\n") ; printf ("Using OGLE light curve format\n") ; #else strcpy (formatStr, "%lf %lf %lf\n") ; // use the default LC format #endif #ifdef DO_PRINT_OmC_PLOT sprintf (filenameSM, "%s.%s.sm", removePath(argv[1]), DO_PRINT_OmC_PLOT) ; if (!(foutSM = fopen (filenameSM, "wt"))) { printf ("ERROR: Couldn't open the output SM script\n") ; return (4) ; } #endif while (10 == fscanf (fin, "%s %lf %lf %*f %lf %*f %lf %*f %lf %*f %lf %*f %lf %*f %lf %*f %lf %*f %*d %*d %*f %*f %*f %*f %*f %*f %*f %*f %*f %*f %*f\n", filename, &period, &e, &r1, &r2, &B1, &B2, &sin_i, &time0, &omega)) { filenameNoPath = removePath(filename) ; B1 = pow(10.0, -0.4 * B1) / integrateWholeDisk(r1) ; B2 = pow(10.0, -0.4 * B2) / integrateWholeDisk(r2) ; omega = omega * M_PI / 180.0 ; printf ("%s (P= %.3f days)\n", filename, period) ; finLC = fopen (filename, "rt") ; if (!finLC) { printf ("Couldn't open LC file\n") ; continue ; } size = 0 ; while (3 == fscanf(finLC, formatStr, &tmpTime, &tmpMag, &tmpErr)) if ((tmpMag < MAX_MAGNITUDE) && (tmpErr > 0.0)) // Sign of an invalid magnitude size++ ; time = (double*)malloc (size * sizeof(double)) ; amp = (double*)malloc (size * sizeof(double)) ; err = (double*)malloc (size * sizeof(double)) ; tmpArr = (double*)malloc (size * sizeof(double)) ; OmCarr = (struct OmCdata*)malloc (size * sizeof(struct OmCdata)) ; if ((!time) || (!amp) || (!err) || (!tmpArr) || (!OmCarr)) { if (time) free(time) ; if (amp) free(amp) ; if (err) free(err) ; if (tmpArr) free(tmpArr) ; if (OmCarr) free(OmCarr) ; printf ("ERROR: Not enough memory") ; fclose(finLC) ; continue ; } //------------------------- rewind(finLC) ; i = 0 ; while ((i < size) && (3 == fscanf(finLC, formatStr, &tmpTime, &tmpMag, &tmpErr))) if ((tmpMag < MAX_MAGNITUDE) && (tmpErr > 0.0)) // Sign of an invalid magnitude { time[i] = tmpTime / period ; amp[i] = pow(10.0, -0.4 * tmpMag) ; // Convert magnitude (logarithmic) to amplitude (linear) err[i] = amp[i] * tmpErr / MAGNITUDE_ERROR_FACTOR ; i++ ; } if (i != size) { printf ("ERROR: Number mismatch\n") ; free(time) ; free(amp) ; free(err) ; fclose (finLC) ; continue ; } i = sortSamples (time, amp, err, size) ; if (i) { printf ("ERROR: Not enough memory for sortSamlpes()\n") ; free(time) ; free(amp) ; free(err) ; fclose (finLC) ; continue ; } //------------------ Find dips --------------------------- tDip1 = getDipCenter (mod1(time0 - (omega / (2.0 * M_PI)) + 1.25), time0, sin_i, e, omega) ; tDip1Wdown = getDipHalfWidth (tDip1, -DIP_FIND_STEP_SIZE, time0, sin_i, e, omega, r1 + r2) ; tDip1Wup = getDipHalfWidth (tDip1, DIP_FIND_STEP_SIZE, time0, sin_i, e, omega, r1 + r2) ; tDip2 = getDipCenter (mod1(time0 - (omega / (2.0 * M_PI)) + 1.75), time0, sin_i, e, omega) ; tDip2Wdown = getDipHalfWidth (tDip2, -DIP_FIND_STEP_SIZE, time0, sin_i, e, omega, r1 + r2) ; tDip2Wup = getDipHalfWidth (tDip2, DIP_FIND_STEP_SIZE, time0, sin_i, e, omega, r1 + r2) ; if (flux (tDip1, time0, sin_i, r1, r2, B1, B2, e, omega, INTEGRATION_STEP) > flux (tDip2, time0, sin_i, r1, r2, B1, B2, e, omega, INTEGRATION_STEP)) { swap (&tDip1, &tDip2) ; swap (&tDip1Wdown, &tDip2Wdown) ; swap (&tDip1Wup, &tDip2Wup) ; } // The plateau lengths should be equal for circular orbits if (tDip2 > tDip1) { lenPlateau12 = ((tDip2 - tDip1) - tDip1Wup - tDip2Wdown) * FRACTION_OF_PLATEAU_TO_USE ; lenPlateau21 = (1.0 - (tDip2 - tDip1) - tDip1Wdown - tDip2Wup) * FRACTION_OF_PLATEAU_TO_USE ; } else // (tDip1 >= tDip2) { lenPlateau12 = (1.0 - (tDip1 - tDip2) - tDip1Wup - tDip2Wdown) * FRACTION_OF_PLATEAU_TO_USE ; lenPlateau21 = ((tDip1 - tDip2) - tDip1Wdown - tDip2Wup) * FRACTION_OF_PLATEAU_TO_USE ; } #ifdef DEBUG printf ("dip1 = %f %f %f (flux=%e)\n", tDip1, tDip1Wup, tDip1Wdown, flux (tDip1, time0, sin_i, r1, r2, B1, B2, e, omega, INTEGRATION_STEP)) ; printf ("dip2 = %f %f %f (flux=%e)\n", tDip2, tDip2Wup, tDip2Wdown, flux (tDip2, time0, sin_i, r1, r2, B1, B2, e, omega, INTEGRATION_STEP)) ; printf ("plateau lengths= %lf %lf\n", lenPlateau12, lenPlateau21) ; #endif //------------------ Find delta-time --------------------------- OmCsize = 0 ; #ifdef DEBUG printf ("dip1:\n") ; #endif for (tDip = tDip1 + floor(time[0]) ; tDip < time[size-1] ; tDip += 1.0) // time[] array is sorted { deltaLocalDip = getLocalDipShift (tDip, tDip1Wdown, tDip1Wup, time, amp, err, tmpArr, size, time0, sin_i, r1, r2, B1, B2, e, omega, lenPlateau21, lenPlateau12, &uncertaintyUp, &uncertaintyDown, &AGCfactor) ; if (deltaLocalDip < MAXFLOAT) { #ifdef DEBUG printf ("%f %f (%f %f) %e\n", period * tDip, OmC_MULTIPLIER * period * deltaLocalDip, OmC_MULTIPLIER * period * uncertaintyUp, OmC_MULTIPLIER * period * uncertaintyDown, -2.5 * log10(AGCfactor)) ; #endif OmCarr[OmCsize].time = tDip ; OmCarr[OmCsize].OmC = deltaLocalDip ; OmCarr[OmCsize].OmC_errUp = uncertaintyUp ; OmCarr[OmCsize].OmC_errDown = uncertaintyDown ; OmCarr[OmCsize].AGCfactor = AGCfactor ; OmCarr[OmCsize].eclipseTag = 1 ; OmCsize++ ; } } #ifdef DEBUG printf ("dip2:\n") ; #endif for (tDip = tDip2 + floor(time[0]) ; tDip < time[size-1] ; tDip += 1.0) // time[] array is sorted { deltaLocalDip = getLocalDipShift (tDip, tDip2Wdown, tDip2Wup, time, amp, err, tmpArr, size, time0, sin_i, r1, r2, B1, B2, e, omega, lenPlateau12, lenPlateau21, &uncertaintyUp, &uncertaintyDown, &AGCfactor) ; if (deltaLocalDip < MAXFLOAT) { #ifdef DEBUG printf ("%f %f (%f %f) %e\n", period * tDip, OmC_MULTIPLIER * period * deltaLocalDip, OmC_MULTIPLIER * period * uncertaintyUp, OmC_MULTIPLIER * period * uncertaintyDown, -2.5 * log10(AGCfactor)) ; #endif OmCarr[OmCsize].time = tDip ; OmCarr[OmCsize].OmC = deltaLocalDip ; OmCarr[OmCsize].OmC_errUp = uncertaintyUp ; OmCarr[OmCsize].OmC_errDown = uncertaintyDown ; OmCarr[OmCsize].AGCfactor = AGCfactor ; OmCarr[OmCsize].eclipseTag = 2 ; OmCsize++ ; } } //---------------- sigma clipping and statistics ------------------------- chi2nu_linear = sigmaClipping (OmCarr, &OmCsize, &slope, &alarmScore) ; if (chi2nu_linear == -1.0) { printf ("ERROR: Too few data points (%d)\n", OmCsize) ; free(time) ; free(amp) ; free(err) ; free (tmpArr) ; free (OmCarr) ; fclose (finLC) ; continue ; } parabFtest = calcParabola_Ftest (filenameNoPath, period, OmCarr, OmCsize, &Pdot, filenameBestFitParab) ; sinusoidalFtest = calcSinusoidal_Ftest (filenameNoPath, period, OmCarr, OmCsize, &bestFreq, &bestFreqAmp, filenameBestFitSin) ; printf ("AlarmScore= %.3lf corrected period= %.12lf N= %d chi2nu_linear = %lf\n", alarmScore, period * (1.0 + slope), OmCsize, chi2nu_linear) ; printf ("parab F-test=%.3lf Pdot= %.3e sec/sec\n", parabFtest, Pdot) ; printf ("sin F-test=%.3lf period= %.3lf days amp= %.3lf sec\n", sinusoidalFtest, 1.0 / bestFreq, OmC_MULTIPLIER * period * bestFreqAmp) ; fprintf (fout, "%s %.12lf %lf %.12lf %d %lf %lf %e %lf %lf %lf\n", filenameNoPath, period, alarmScore, period * (1.0 + slope), OmCsize, chi2nu_linear, parabFtest, Pdot, sinusoidalFtest, 1.0 / bestFreq, OmC_MULTIPLIER * period * bestFreqAmp) ; fflush (fout) ; #ifdef DO_PRINT_OmC_PLOT hasPrimaryEclipses = 0 ; for (i = 0 ; (i < OmCsize) && (!hasPrimaryEclipses) ; i++) if (OmCarr[i].eclipseTag == 1) hasPrimaryEclipses = 1 ; hasSecondaryEclipses = 0 ; for (i = 0 ; (i < OmCsize) && (!hasSecondaryEclipses) ; i++) if (OmCarr[i].eclipseTag == 2) hasSecondaryEclipses = 1 ; if (hasPrimaryEclipses) { sprintf (filenamePlot1, "%s_P%.3f.%s.1", filenameNoPath, period, DO_PRINT_OmC_PLOT) ; foutPlot1 = fopen (filenamePlot1, "wt") ; if (!foutPlot1) printf ("ERROR: Couldn't open one of the output plot files\n") ; else { for (i = 0 ; i < OmCsize ; i++) if (OmCarr[i].eclipseTag == 1) fprintf (foutPlot1, "%lf %lf %lf %lf %e\n", period * OmCarr[i].time, OmC_MULTIPLIER * period * OmCarr[i].OmC, OmC_MULTIPLIER * period * OmCarr[i].OmC_errUp, OmC_MULTIPLIER * period * OmCarr[i].OmC_errDown, -2.5 * log10(OmCarr[i].AGCfactor)) ; fclose (foutPlot1) ; } } if (hasSecondaryEclipses) { sprintf (filenamePlot2, "%s_P%.3f.%s.2", filenameNoPath, period, DO_PRINT_OmC_PLOT) ; foutPlot2 = fopen (filenamePlot2, "wt") ; if (!foutPlot2) printf ("ERROR: Couldn't open one of the output plot files\n") ; else { for (i = 0 ; i < OmCsize ; i++) if (OmCarr[i].eclipseTag == 2) fprintf (foutPlot2, "%lf %lf %lf %lf %e\n", period * OmCarr[i].time, OmC_MULTIPLIER * period * OmCarr[i].OmC, OmC_MULTIPLIER * period * OmCarr[i].OmC_errUp, OmC_MULTIPLIER * period * OmCarr[i].OmC_errDown, -2.5 * log10(OmCarr[i].AGCfactor)) ; fclose (foutPlot2) ; } } // for the readability of the O-C plot's X-axis baseTime = (long)(floor(OmCarr[0].time * period / BASE_TIME_ROUNDING_FACTOR) * BASE_TIME_ROUNDING_FACTOR) ; writeSM (foutSM, filenameSM, filenameNoPath, filenamePlot1, filenamePlot2, period, baseTime, hasPrimaryEclipses, hasSecondaryEclipses, filenameBestFitParab, filenameBestFitSin) ; #endif free(time) ; free(amp) ; free(err) ; free (tmpArr) ; free (OmCarr) ; fclose (finLC) ; } fclose (fin) ; fclose (fout) ; #ifdef DO_PRINT_OmC_PLOT fprintf (foutSM, "hardcopy\n") ; fclose (foutSM) ; #endif return (0) ; } /* Todo: 1. better y-axis tick labels (>= 10,000 ; <= -10,000) ??? 2. remove long-term trends (in DEBiL, iterative) ? */