/********************************************************** * * This program simulates light curves, for 31-parameter (raw - direct from DEBiL program) * * compile: gcc timing2periodCorrect.c -o timing2periodCorrect -Wall -lm -O4 * * run: timing2periodCorrect cmDra2.dbl cmDra2.timing cmDra2.timing_err 0.4465 0.1855 cmDra2.period cmDra2.period_ecc * **********************************************************/ #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_BEST_FIT_LINEAR "linear" #define DO_PRINT_BEST_FIT_OFFSET1 "offset1" #define DO_PRINT_BEST_FIT_OFFSET2 "offset2" // The thresholds used to determine if the EB has an eccentric orbit #define THRESHOLD_F_TEST 0.5 // the larger the F-test value is, the more systems will be considered eccentric #define THRESHOLD_E_COSW 0.002 //#define RESET_ECC_TO_ZERO // if there are many systems with reported eccentrisities that you do not believe //#define DEBUG #define OUTLIER_SIGMA_LIMIT 6.0 // for "sigma clipping" #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 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 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 *filenameBestFitLinear, char *filenameBestFitOffset1, char *filenameBestFitOffset2) { 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 (*filenameBestFitLinear) { fprintf (foutSM, "ltype 1\n") ; fprintf (foutSM, "ctype blue\n") ; fprintf (foutSM, "data %s\n", filenameBestFitLinear) ; 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 (*filenameBestFitOffset1) { fprintf (foutSM, "ltype 2\n") ; fprintf (foutSM, "ctype red\n") ; fprintf (foutSM, "data %s\n", filenameBestFitOffset1) ; 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 (*filenameBestFitOffset2) { fprintf (foutSM, "ltype 2\n") ; fprintf (foutSM, "ctype green\n") ; fprintf (foutSM, "data %s\n", filenameBestFitOffset2) ; 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 (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 ------------------------ // Removes the outliers (from a weighted first-order regression) from the input data array. // Finds the best linear fit, and optionally plots it. // NOTE: The O-C curve is NOT normalized (de-trended) // Returns the final reduced chi-squared (chi2nu) around the linear regression // Upon error returns -1.0 (too few data points, or infinite slope) // Note that if there is an error, the values of (*pslope) is set to zero double sigmaClipping (struct OmCdata *OmCarr, int *pOmCsize, double *pslope, double period, char *filenameNoPath, char *filenameBestFitLinear) { int i, doClip, numClipped = 0, OmCsize = *pOmCsize ; double N, A, B, C, D, E, tmp, slope, offset, chi2nu, t, x, invSigma2, tAvr ; double sigma2, diff ; FILE *foutLinear ; *pslope = 0.0 ; do { if (OmCsize < 2) { 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)) ; if (OmCsize > 2) chi2nu /= (OmCsize-2) ; else chi2nu = MAXFLOAT ; #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 ; qsort (OmCarr, OmCsize, sizeof(struct OmCdata), comparOmCdata) ; // ------------- plot best linear fit ---------------------- #ifdef DO_PRINT_BEST_FIT_LINEAR sprintf (filenameBestFitLinear, "%s_P%.3f.%s", filenameNoPath, period, DO_PRINT_BEST_FIT_LINEAR) ; foutLinear = fopen (filenameBestFitLinear, "wt") ; if (!foutLinear) { printf ("ERROR: Couln't open '%s' for writing best linear fit\n", filenameBestFitLinear) ; *filenameBestFitLinear = 0 ; } else { fprintf (foutLinear, "%lf %lf\n", OmCarr[0].time * period, OmC_MULTIPLIER * period * ((OmCarr[0].time * slope) + offset)) ; fprintf (foutLinear, "%lf %lf\n", OmCarr[OmCsize-1].time * period, OmC_MULTIPLIER * period * ((OmCarr[OmCsize-1].time * slope) + offset)) ; fclose (foutLinear) ; } #else *filenameBestFitLinear = 0 ; #endif return (chi2nu) ; } // Finds the best linear fit for the primary/secondary eclipses (same slope), and optionally plots them. // Using this result is estimates the value of e*cos(omega), where e=eccentricity and omega=angle of perihelion // Returns the final reduced chi-squared (chi2nu) around the two linear regressions // Upon error returns -1.0 (too few data points, or infinite slope) or MAXFLOAT (N=3) data points // Note that if there is an error, the values of (*pslope) and (*p_ecosw) are not changed double calc_ecosw (struct OmCdata *OmCarr, int OmCsize, double *pslope, double period, double *p_ecosw, char *filenameNoPath, char *filenameBestFitOffset1, char *filenameBestFitOffset2) { int i ; double N1, A1, B1, C1, D1, E1, N2, A2, B2, C2, D2, E2 ; double tmp, slope, offset1, offset2, chi2nu, t, x, invSigma2, tAvr ; FILE *foutLinear ; if (OmCsize < 3) return (-1.0) ; N1 = 0.0 ; A1 = 0.0 ; B1 = 0.0 ; C1 = 0.0 ; D1 = 0.0 ; E1 = 0.0 ; N2 = 0.0 ; A2 = 0.0 ; B2 = 0.0 ; C2 = 0.0 ; D2 = 0.0 ; E2 = 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 if (OmCarr[i].eclipseTag == 1) { N1 += invSigma2 ; A1 += (t * invSigma2) ; B1 += (t * t * invSigma2) ; C1 += (t * x * invSigma2) ; D1 += (x * invSigma2) ; E1 += (x * x * invSigma2) ; } else // (OmCarr[i].eclipseTag == 2) { N2 += invSigma2 ; A2 += (t * invSigma2) ; B2 += (t * t * invSigma2) ; C2 += (t * x * invSigma2) ; D2 += (x * invSigma2) ; E2 += (x * x * invSigma2) ; } } if ((N1 == 0.0) || (N2 == 0.0)) return (-1.0) ; tmp = (N1*A2*A2) + (N2*A1*A1) - (N1*N2*(B1+B2)) ; if ((tmp < EPSILON) && (tmp > -EPSILON)) return (-1.0) ; slope = ((N1*A2*D2) + (N2*A1*D1) - (N1*N2*(C1+C2))) / tmp ; offset1 = ((D1*A2*A2) - (D2*A1*A2) - (D1*N2*(B1+B2)) + (A1*N2*(C1+C2))) / tmp ; offset2 = ((D2*A1*A1) - (D1*A1*A2) - (D2*N1*(B1+B2)) + (A2*N1*(C1+C2))) / tmp ; chi2nu = (slope*slope*(B1+B2)) + (N1*offset1*offset1) + (N2*offset2*offset2) + E1 + E2 + (2.0*slope*((offset1*A1) + (offset2*A2) - C1 - C2)) - (2.0*offset1*D1) - (2.0*offset2*D2) ; if (OmCsize > 3) chi2nu /= (OmCsize-3) ; else chi2nu = MAXFLOAT ; *pslope = slope ; *p_ecosw += (0.5 * M_PI * (offset2 - offset1)) ; // no need to divide by period, since offset1,2 are already nomalized offset1 -= (slope * tAvr) ; offset2 -= (slope * tAvr) ; #ifdef DEBUG printf ("Ecc Num= %d (N=%e %e) slope= %lf offset1= %lf offset2= %lf e*cos(w)= %lf delta e*cos(w)= %lf\n", OmCsize, N1, N2, slope, offset1, offset2, *p_ecosw, 0.5 * M_PI * (offset2 - offset1)) ; #endif // ------------- plot best linear fit ---------------------- #ifdef DO_PRINT_BEST_FIT_OFFSET1 sprintf (filenameBestFitOffset1, "%s_P%.3f.%s", filenameNoPath, period, DO_PRINT_BEST_FIT_OFFSET1) ; foutLinear = fopen (filenameBestFitOffset1, "wt") ; if (!foutLinear) { printf ("ERROR: Couln't open '%s' for writing best linear-offset 1 fit\n", filenameBestFitOffset1) ; *filenameBestFitOffset1 = 0 ; } else { fprintf (foutLinear, "%lf %lf\n", OmCarr[0].time * period, OmC_MULTIPLIER * period * ((OmCarr[0].time * slope) + offset1)) ; fprintf (foutLinear, "%lf %lf\n", OmCarr[OmCsize-1].time * period, OmC_MULTIPLIER * period * ((OmCarr[OmCsize-1].time * slope) + offset1)) ; fclose (foutLinear) ; } #else *filenameBestFitOffset1 = 0 ; #endif #ifdef DO_PRINT_BEST_FIT_OFFSET2 sprintf (filenameBestFitOffset2, "%s_P%.3f.%s", filenameNoPath, period, DO_PRINT_BEST_FIT_OFFSET2) ; foutLinear = fopen (filenameBestFitOffset2, "wt") ; if (!foutLinear) { printf ("ERROR: Couln't open '%s' for writing best linear-offset 2 fit\n", filenameBestFitOffset2) ; *filenameBestFitOffset2 = 0 ; } else { fprintf (foutLinear, "%lf %lf\n", OmCarr[0].time * period, OmC_MULTIPLIER * period * ((OmCarr[0].time * slope) + offset2)) ; fprintf (foutLinear, "%lf %lf\n", OmCarr[OmCsize-1].time * period, OmC_MULTIPLIER * period * ((OmCarr[OmCsize-1].time * slope) + offset2)) ; fclose (foutLinear) ; } #else *filenameBestFitOffset2 = 0 ; #endif return (chi2nu) ; } int main (int argc, char **argv) { FILE *fin, *finLC, *fout, *ferr, *fPeriod = 0, *fPeriodEcc = 0, *fUncorrectedDebil = 0 ; 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, new_ecc ; double x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18 ; double x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31 ; int x19, x20 ; char filename[MAX_STRING_LEN], *filenameNoPath ; char filenameBestFitLinear[MAX_STRING_LEN] ; char filenameBestFitOffset1[MAX_STRING_LEN], filenameBestFitOffset2[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, chi2nu_linear, chi2nu_offsets ; double uncertaintyUp, uncertaintyDown, AGCfactor, slope, ecosw, ecosw_orig ; 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 != 4) && (argc != 6) && (argc != 8) && (argc != 9)) { printf ("Usage: %s [ ] [<(non-ecc EB period list> [uncorrected DEBiL systems]]\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") ; fclose (fin) ; return (3) ; } if (!(ferr = fopen (argv[3], "wt"))) { printf ("ERROR: Couldn't open the output error file\n") ; fclose (fin) ; fclose (fout) ; return (4) ; } if (argc >= 5) { quadA = atof(argv[4]) ; quadB = atof(argv[5]) ; } printf ("Limb darkening coefficients: A=%f ; B=%f\n", quadA, quadB) ; setLimbDarkening (quadA, quadB) ; if (argc >= 7) { if (!(fPeriod = fopen (argv[6], "wt"))) { printf ("ERROR: Couldn't open the output non-ecc period file\n") ; fclose (fin) ; fclose (fout) ; fclose (ferr) ; return (5) ; } if (!(fPeriodEcc = fopen (argv[7], "wt"))) { printf ("ERROR: Couldn't open the output non-ecc period file\n") ; fclose (fin) ; fclose (fout) ; fclose (ferr) ; fclose (fPeriod) ; return (6) ; } } if (argc >= 9) { if (!(fUncorrectedDebil = fopen (argv[8], "wt"))) { printf ("ERROR: Couldn't open the output uncorrected DEBiL file\n") ; fclose (fin) ; fclose (fout) ; fclose (ferr) ; fclose (fPeriod) ; fclose (fPeriodEcc) ; return (7) ; } } #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") ; fclose (fin) ; fclose (fout) ; fclose (ferr) ; fclose (fPeriod) ; fclose (fPeriodEcc) ; return (8) ; } #endif while (31 == fscanf (fin, "%s %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", filename, &x2, &x3, &x4, &x5, &x6, &x7, &x8, &x9, &x10, &x11, &x12, &x13, &x14, &x15, &x16, &x17, &x18, &x19, &x20, &x21, &x22, &x23, &x24, &x25, &x26, &x27, &x28, &x29, &x30, &x31)) { filenameNoPath = removePath(filename) ; period = x2 ; #ifdef RESET_ECC_TO_ZERO e = 0.0 ; #else e = x3 ; #endif r1 = x5 ; r2 = x7 ; B1 = pow(10.0, -0.4 * x9) / integrateWholeDisk(r1) ; B2 = pow(10.0, -0.4 * x11) / integrateWholeDisk(r2) ; sin_i = x13 ; time0 = x15 ; omega = x17 * 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") ; fprintf (ferr, "%s %.12lf [Couldn't open LC file]\n", filename, period) ; if (fUncorrectedDebil) fprintf (fUncorrectedDebil, "%s %.12f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %d %d %f %f %f %f %f %f %f %f %f %f %f\n", filename, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31) ; 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") ; fprintf (ferr, "%s %.12lf [Not enough memory]\n", filename, period) ; if (fUncorrectedDebil) fprintf (fUncorrectedDebil, "%s %.12f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %d %d %f %f %f %f %f %f %f %f %f %f %f\n", filename, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31) ; 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") ; fprintf (ferr, "%s %.12lf [Number mismatch]\n", filename, period) ; if (fUncorrectedDebil) fprintf (fUncorrectedDebil, "%s %.12f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %d %d %f %f %f %f %f %f %f %f %f %f %f\n", filename, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31) ; 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") ; fprintf (ferr, "%s %.12lf [Not enough memory for sortSamlpes()]\n", filename, period) ; if (fUncorrectedDebil) fprintf (fUncorrectedDebil, "%s %.12f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %d %d %f %f %f %f %f %f %f %f %f %f %f\n", filename, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31) ; 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, period, filenameNoPath, filenameBestFitLinear) ; if (chi2nu_linear == -1.0) { printf ("ERROR: Too few data points for sigmaClipping() [N=%d]\n", OmCsize) ; fprintf (ferr, "%s %.12lf [Too few data points (%d)]\n", filename, period, OmCsize) ; if (fUncorrectedDebil) fprintf (fUncorrectedDebil, "%s %.12f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %d %d %f %f %f %f %f %f %f %f %f %f %f\n", filename, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31) ; free(time) ; free(amp) ; free(err) ; free (tmpArr) ; free (OmCarr) ; fclose (finLC) ; continue ; } #ifdef DEBUG printf ("linear slope= %lf\n", slope) ; #endif ecosw = ecosw_orig = e * cos(omega) ; chi2nu_offsets = calc_ecosw (OmCarr, OmCsize, &slope, period, &ecosw, filenameNoPath, filenameBestFitOffset1, filenameBestFitOffset2) ; if (chi2nu_offsets == -1.0) { printf ("Warning: Too few data points for calc_ecosw() [N=%d]\n", OmCsize) ; *filenameBestFitOffset1 = 0 ; *filenameBestFitOffset2 = 0 ; } if (e == 0.0) new_ecc = fabs(ecosw) ; else new_ecc = ecosw / cos(omega) ; //-------------------------- printf ("Orig period= %.12lf corrected period= %.12lf e*cos(w)= %lf [e= %lf] N= %d ", period, period * (1.0 + slope), ecosw, new_ecc, OmCsize) ; if ((chi2nu_linear >= 0.0) && (chi2nu_linear < MAXFLOAT)) printf ("chi2nu_linear= %lf ", chi2nu_linear) ; else printf ("chi2nu_linear= N/A ") ; if ((chi2nu_offsets >= 0.0) && (chi2nu_offsets < MAXFLOAT)) printf ("chi2nu_offsets= %lf\n", chi2nu_offsets) ; else printf ("chi2nu_offsets= N/A\n") ; //---- fprintf (fout, "%s %.12lf %.12lf %lf %lf %d %lf ", filename, period, period * (1.0 + slope), ecosw, new_ecc, OmCsize, chi2nu_linear) ; if ((chi2nu_linear >= 0.0) && (chi2nu_linear < MAXFLOAT)) fprintf (fout, "%lf ", chi2nu_linear) ; else fprintf (fout, "-1 ") ; if ((chi2nu_offsets >= 0.0) && (chi2nu_offsets < MAXFLOAT)) fprintf (fout, "%lf\n", chi2nu_offsets) ; else fprintf (fout, "-1\n") ; fflush (fout) ; //---- if (fPeriod) { if ((fabs(ecosw) < THRESHOLD_E_COSW) || ((fabs(ecosw_orig) < THRESHOLD_E_COSW) && ((chi2nu_offsets == -1.0) || (chi2nu_offsets > (chi2nu_linear * THRESHOLD_F_TEST))))) { fprintf (fPeriod, "%s %.12lf\n", filename, period * (1.0 + slope)) ; fflush (fPeriod) ; } else { fprintf (fPeriodEcc, "%s %.12lf\n", filename, period * (1.0 + slope)) ; fflush (fPeriodEcc) ; } } #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, filenameBestFitLinear, filenameBestFitOffset1, filenameBestFitOffset2) ; #endif free(time) ; free(amp) ; free(err) ; free (tmpArr) ; free (OmCarr) ; fclose (finLC) ; } #ifdef DO_PRINT_OmC_PLOT fprintf (foutSM, "hardcopy\n") ; fclose (foutSM) ; printf ("\nCreated files:\n") ; printf (" *.plot.1 - O-C data points from the primary eclipses\n") ; printf (" *.plot.2 - O-C data points from the secondary eclipses\n") ; printf (" *.linear - A linear fit to all data points\n") ; printf (" *.offset1 - An offsetted linear fit to the primary eclipse data points\n") ; printf (" *.offset2 - An offsetted linear fit to the secondary eclipse data points\n") ; printf (" %s - The SM file (a SuperMongo script)\n", filenameSM) ; printf ("\nTo run the SM file, enter: 'sm', then 'inp %s', then 'quit'.\n", filenameSM) ; printf ("This will create a postscript file ('%s.ps') of the model fits catalogue.\n\n", filenameSM) ; #endif fclose (fin) ; fclose (fout) ; fclose (ferr) ; if (fPeriod) fclose (fPeriod) ; if (fPeriodEcc) fclose (fPeriodEcc) ; if (fUncorrectedDebil) fclose (fUncorrectedDebil) ; return (0) ; } /* Todo: 1. better y-axis tick labels (>= 10,000 ; <= -10,000) ??? 2. remove long-term trends (in DEBiL, iterative) ? */