diff --git a/smoke_test/assets/sources/smoke_source_binary_catalog.dat b/smoke_test/assets/sources/smoke_source_binary_catalog.dat index 8ed3b04..df08dd9 100644 --- a/smoke_test/assets/sources/smoke_source_binary_catalog.dat +++ b/smoke_test/assets/sources/smoke_source_binary_catalog.dat @@ -1,7 +1,7 @@ -R062 Z087 Y106 J129 W146 H158 F184 K213 2MASS_Ks 2MASS_J 2MASS_H Bessell_U Bessell_B Bessell_V Bessell_I Bessell_R Kepler_Kp TESS DECam_z DECam_u DECam_g DECam_r DECam_i DECam_Y Gaia_G_EDR3 Gaia_BP_EDR3 Gaia_RP_EDR3 VISTA_Z VISTA_Y VISTA_J VISTA_H VISTA_Ks mul mub Vr U V W iMass CL age Teff logg pop Mass Mbol Radius [Fe/H] l b RA2000.0 DEC2000.0 Dist x y z A_Ks [alpha/Fe] ID Is_Binary primary_ID combined_logP -2.5212534e+01 2.3437151e+01 2.2946926e+01 2.2698885e+01 2.2737316e+01 2.2564209e+01 2.2718139e+01 2.2750813e+01 2.0916308e+01 2.1813254e+01 2.1137727e+01 2.8960719e+01 2.7692111e+01 2.6146455e+01 2.3327938e+01 2.4879330e+01 2.4833463e+01 2.3528349e+01 2.3269717e+01 2.8460170e+01 2.6830101e+01 2.5344168e+01 2.3886079e+01 2.3095472e+01 2.4810155e+01 2.6388821e+01 2.3603958e+01 2.2911913e+01 2.2447639e+01 2.1762028e+01 2.1155567e+01 2.0943504e+01 -1.0430274e+01 2.9710474e+00 5.3098679e+01 6.6088945e+01 1.0138185e+02 4.8754663e+01 1.3530942e-01 0.0000000e+00 1.0000000e+01 3.1466277e+03 5.1644453e+00 0.0000000e+00 1.3530814e-01 1.1389806e+01 1.5813397e-01 -9.4970712e-02 3.2006745e-03 4.2958514e-03 2.6640150e+02 -2.8931190e+01 2.9168337e+00 -5.2611722e+00 1.6294106e-04 1.1155346e-02 0.0000000e+00 0 1 1 1 1.0 -1.8280023e+01 1.8068850e+01 1.8049266e+01 1.8055254e+01 1.8117239e+01 1.8116642e+01 1.8302014e+01 1.8562267e+01 1.6727762e+01 1.7164161e+01 1.6776748e+01 1.9448572e+01 1.9185598e+01 1.8437426e+01 1.7656769e+01 1.8021182e+01 1.8199912e+01 1.7754535e+01 1.8063731e+01 2.0042738e+01 1.8749093e+01 1.8215518e+01 1.8096851e+01 1.8045478e+01 1.8284698e+01 1.8650127e+01 1.7755000e+01 1.7576242e+01 1.7436794e+01 1.7135492e+01 1.6781831e+01 1.6731542e+01 5.4187211e-02 3.8186918e+00 -1.1389063e+02 -1.0081186e+02 2.4678124e+02 9.1264341e+01 9.0224882e-01 0.0000000e+00 1.0000000e+01 5.4611661e+03 4.4174184e+00 0.0000000e+00 9.0202770e-01 5.0507352e+00 9.7265615e-01 6.4701074e-02 3.5999997e+02 -3.7644077e-03 2.6640744e+02 -2.8938149e+01 4.5983030e+00 -3.5797075e+00 -2.6445424e-06 7.1391942e-03 0.0000000e+00 0 2 2 1 1.0 -2.7189580e+01 2.4975054e+01 2.4260565e+01 2.3968756e+01 2.4030988e+01 2.3845912e+01 2.4039567e+01 2.3989889e+01 2.2155384e+01 2.3048865e+01 2.2390443e+01 3.1118061e+01 3.0030954e+01 2.8490833e+01 2.5036992e+01 2.6885617e+01 2.6681963e+01 2.5118838e+01 2.4731204e+01 3.0059395e+01 2.9036582e+01 2.7540074e+01 2.5635858e+01 2.4463607e+01 2.6583638e+01 2.8723106e+01 2.5259802e+01 2.4469663e+01 2.3817850e+01 2.2989283e+01 2.2414511e+01 2.2190391e+01 -6.5516244e+00 1.1716060e+00 -7.9343188e+01 -6.6373337e+01 9.2833791e+01 3.5265923e+01 1.4345340e-01 0.0000000e+00 1.0000000e+01 2.9315876e+03 5.1305404e+00 0.0000000e+00 1.4345205e-01 1.1600987e+01 1.6530046e-01 1.6760650e-01 4.3851388e-03 -2.5059108e-03 2.6640884e+02 -2.8933722e+01 4.9185785e+00 -3.2594326e+00 3.7644395e-04 6.5604172e-03 0.0000000e+00 0 3 0 0 0.0 -2.6432685e+01 2.4491108e+01 2.3915671e+01 2.3656481e+01 2.3706363e+01 2.3533906e+01 2.3701900e+01 2.3710237e+01 2.1875732e+01 2.2758205e+01 2.2099233e+01 3.0094746e+01 2.8982746e+01 2.7485130e+01 2.4443064e+01 2.6113111e+01 2.6006401e+01 2.4604554e+01 2.4293562e+01 2.9474944e+01 2.8101540e+01 2.6642829e+01 2.5017614e+01 2.4081378e+01 2.5955676e+01 2.7725745e+01 2.4703953e+01 2.3972942e+01 2.3432310e+01 2.2704340e+01 2.2119082e+01 2.1905447e+01 -7.7780939e+00 -2.3327942e-01 -2.3962213e+01 -1.1064623e+01 6.4507357e+01 2.3963224e+00 1.5178442e-01 0.0000000e+00 1.0000000e+01 3.0756073e+03 5.1140763e+00 0.0000000e+00 1.5178288e-01 1.1246362e+01 1.7682473e-01 8.4404856e-02 2.6359830e-03 5.2891509e-03 2.6640020e+02 -2.8931154e+01 4.9113800e+00 -3.2666297e+00 2.2595581e-04 7.2438851e-03 0.0000000e+00 0 4 0 0 0.0 -4.7021105e+01 3.4017351e+01 2.8549183e+01 2.4842273e+01 2.4605516e+01 2.2513490e+01 2.1428106e+01 2.0791151e+01 1.8956646e+01 2.4518952e+01 2.0729039e+01 7.5145964e+01 6.6431818e+01 5.2662162e+01 3.7022015e+01 4.5436628e+01 4.8035999e+01 4.0072213e+01 3.2192061e+01 7.2625737e+01 6.1205873e+01 4.5835721e+01 3.8642084e+01 2.9947132e+01 4.4114127e+01 5.4168672e+01 3.7180335e+01 3.2877854e+01 2.8552558e+01 2.4147819e+01 2.0674371e+01 1.8947702e+01 -4.6698104e+00 -2.4039011e+00 -1.4320101e+02 -1.3046019e+02 1.3330771e+02 -4.9805759e+01 8.1309438e-01 0.0000000e+00 1.0000000e+01 6.0603449e+03 4.3935837e+00 0.0000000e+00 8.1280571e-01 4.6548136e+00 9.4780662e-01 -7.9940091e-01 3.5997144e+02 1.6664060e-02 2.6637053e+02 -2.8951850e+01 5.0758087e+00 -3.1022000e+00 -2.5300458e-03 7.9249568e-03 1.9693590e+01 0 5 1 5 2.0 -4.7021105e+01 3.4017351e+01 2.8549183e+01 2.4842273e+01 2.4605516e+01 2.2513490e+01 2.1428106e+01 2.0791151e+01 1.8956646e+01 2.4518952e+01 2.0729039e+01 7.5145964e+01 6.6431818e+01 5.2662162e+01 3.7022015e+01 4.5436628e+01 4.8035999e+01 4.0072213e+01 3.2192061e+01 7.2625737e+01 6.1205873e+01 4.5835721e+01 3.8642084e+01 2.9947132e+01 4.4114127e+01 5.4168672e+01 3.7180335e+01 3.2877854e+01 2.8552558e+01 2.4147819e+01 2.0674371e+01 1.8947702e+01 -4.6698104e+00 -2.4039011e+00 -1.4320101e+02 -1.3046019e+02 1.3330771e+02 -4.9805759e+01 8.1309438e-01 0.0000000e+00 1.0000000e+01 6.0603449e+03 4.3935837e+00 0.0000000e+00 8.1280571e-01 4.6548136e+00 9.4780662e-01 -7.9940091e-01 3.5997144e+02 1.6664060e-02 2.6637053e+02 -2.8951850e+01 5.0758087e+00 -3.1022000e+00 -2.5300458e-03 7.9249568e-03 1.9693590e+01 0 6 2 5 2.0 +R062 Z087 Y106 J129 W146 H158 F184 K213 2MASS_Ks 2MASS_J 2MASS_H Bessell_U Bessell_B Bessell_V Bessell_I Bessell_R Kepler_Kp TESS DECam_z DECam_u DECam_g DECam_r DECam_i DECam_Y Gaia_G_EDR3 Gaia_BP_EDR3 Gaia_RP_EDR3 VISTA_Z VISTA_Y VISTA_J VISTA_H VISTA_Ks mul mub Vr U V W iMass CL age Teff logg pop Mass Mbol Radius [Fe/H] l b RA2000.0 DEC2000.0 Dist x y z A_Ks [alpha/Fe] ID Is_Binary primary_ID combined_logP +2.5212534e+01 2.3437151e+01 2.2946926e+01 2.2698885e+01 2.2737316e+01 2.2564209e+01 2.2718139e+01 2.2750813e+01 2.0916308e+01 2.1813254e+01 2.1137727e+01 2.8960719e+01 2.7692111e+01 2.6146455e+01 2.3327938e+01 2.4879330e+01 2.4833463e+01 2.3528349e+01 2.3269717e+01 2.8460170e+01 2.6830101e+01 2.5344168e+01 2.3886079e+01 2.3095472e+01 2.4810155e+01 2.6388821e+01 2.3603958e+01 2.2911913e+01 2.2447639e+01 2.1762028e+01 2.1155567e+01 2.0943504e+01 5.4187211e-02 3.8186918e+00 -1.1389063e+02 -1.0081186e+02 2.4678124e+02 9.1264341e+01 1.3530942e-01 0.0000000e+00 1.0000000e+01 3.1466277e+03 5.1644453e+00 0.0000000e+00 1.3530814e-01 1.1389806e+01 1.5813397e-01 -9.4970712e-02 3.2006745e-03 -3.7644077e-03 2.6640150e+02 -2.8931190e+01 2.9168337e+00 -5.2611722e+00 1.6294106e-04 1.1155346e-02 0.0000000e+00 0 1 1 1 1.0 +1.8280023e+01 1.8068850e+01 1.8049266e+01 1.8055254e+01 1.8117239e+01 1.8116642e+01 1.8302014e+01 1.8562267e+01 1.6727762e+01 1.7164161e+01 1.6776748e+01 1.9448572e+01 1.9185598e+01 1.8437426e+01 1.7656769e+01 1.8021182e+01 1.8199912e+01 1.7754535e+01 1.8063731e+01 2.0042738e+01 1.8749093e+01 1.8215518e+01 1.8096851e+01 1.8045478e+01 1.8284698e+01 1.8650127e+01 1.7755000e+01 1.7576242e+01 1.7436794e+01 1.7135492e+01 1.6781831e+01 1.6731542e+01 5.4187211e-02 3.8186918e+00 -1.1389063e+02 -1.0081186e+02 2.4678124e+02 9.1264341e+01 9.0224882e-01 0.0000000e+00 1.0000000e+01 5.4611661e+03 4.4174184e+00 0.0000000e+00 9.0202770e-01 5.0507352e+00 9.7265615e-01 6.4701074e-02 3.5999997e+02 -3.7644077e-03 2.6640744e+02 -2.8938149e+01 4.5983030e+00 -3.5797075e+00 -2.6445424e-06 7.1391942e-03 0.0000000e+00 0 2 2 1 1.0 +2.7189580e+01 2.4975054e+01 2.4260565e+01 2.3968756e+01 2.4030988e+01 2.3845912e+01 2.4039567e+01 2.3989889e+01 2.2155384e+01 2.3048865e+01 2.2390443e+01 3.1118061e+01 3.0030954e+01 2.8490833e+01 2.5036992e+01 2.6885617e+01 2.6681963e+01 2.5118838e+01 2.4731204e+01 3.0059395e+01 2.9036582e+01 2.7540074e+01 2.5635858e+01 2.4463607e+01 2.6583638e+01 2.8723106e+01 2.5259802e+01 2.4469663e+01 2.3817850e+01 2.2989283e+01 2.2414511e+01 2.2190391e+01 -6.5516244e+00 1.1716060e+00 -7.9343188e+01 -6.6373337e+01 9.2833791e+01 3.5265923e+01 1.4345340e-01 0.0000000e+00 1.0000000e+01 2.9315876e+03 5.1305404e+00 0.0000000e+00 1.4345205e-01 1.1600987e+01 1.6530046e-01 1.6760650e-01 4.3851388e-03 -2.5059108e-03 2.6640884e+02 -2.8933722e+01 4.9185785e+00 -3.2594326e+00 3.7644395e-04 6.5604172e-03 0.0000000e+00 0 3 0 0 0.0 +2.6432685e+01 2.4491108e+01 2.3915671e+01 2.3656481e+01 2.3706363e+01 2.3533906e+01 2.3701900e+01 2.3710237e+01 2.1875732e+01 2.2758205e+01 2.2099233e+01 3.0094746e+01 2.8982746e+01 2.7485130e+01 2.4443064e+01 2.6113111e+01 2.6006401e+01 2.4604554e+01 2.4293562e+01 2.9474944e+01 2.8101540e+01 2.6642829e+01 2.5017614e+01 2.4081378e+01 2.5955676e+01 2.7725745e+01 2.4703953e+01 2.3972942e+01 2.3432310e+01 2.2704340e+01 2.2119082e+01 2.1905447e+01 -7.7780939e+00 -2.3327942e-01 -2.3962213e+01 -1.1064623e+01 6.4507357e+01 2.3963224e+00 1.5178442e-01 0.0000000e+00 1.0000000e+01 3.0756073e+03 5.1140763e+00 0.0000000e+00 1.5178288e-01 1.1246362e+01 1.7682473e-01 8.4404856e-02 2.6359830e-03 5.2891509e-03 2.6640020e+02 -2.8931154e+01 4.9113800e+00 -3.2666297e+00 2.2595581e-04 7.2438851e-03 0.0000000e+00 0 4 0 0 0.0 +4.7021105e+01 3.4017351e+01 2.8549183e+01 2.4842273e+01 2.4605516e+01 2.2513490e+01 2.1428106e+01 2.0791151e+01 1.8956646e+01 2.4518952e+01 2.0729039e+01 7.5145964e+01 6.6431818e+01 5.2662162e+01 3.7022015e+01 4.5436628e+01 4.8035999e+01 4.0072213e+01 3.2192061e+01 7.2625737e+01 6.1205873e+01 4.5835721e+01 3.8642084e+01 2.9947132e+01 4.4114127e+01 5.4168672e+01 3.7180335e+01 3.2877854e+01 2.8552558e+01 2.4147819e+01 2.0674371e+01 1.8947702e+01 -9.0335610e+00 -2.7360242e+00 -2.0160392e+02 -1.9107442e+02 5.1533963e+01 -4.5152342e+01 8.1309438e-01 0.0000000e+00 1.0000000e+01 6.0603449e+03 4.3935837e+00 0.0000000e+00 8.1280571e-01 4.6548136e+00 9.4780662e-01 -7.9940091e-01 3.5997144e+02 1.6664060e-02 2.6637053e+02 -2.8951850e+01 5.0758087e+00 -3.1022000e+00 -2.5300458e-03 7.9249568e-03 0.0000000e+00 0 5 1 5 3.0 +3.1793245e+01 2.7449580e+01 2.5641867e+01 2.4639045e+01 2.4657181e+01 2.4083239e+01 2.4042948e+01 2.3798342e+01 2.1963837e+01 2.3810243e+01 2.2533581e+01 0.0000000e+00 0.0000000e+00 3.4035114e+01 2.8110429e+01 3.1294101e+01 3.1417767e+01 2.8566321e+01 2.6860705e+01 3.7798231e+01 3.5635318e+01 3.2055912e+01 2.8903792e+01 2.6126297e+01 3.0666724e+01 3.4440840e+01 2.8305075e+01 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 -9.0335610e+00 -2.7360242e+00 -2.0160392e+02 -1.9107442e+02 5.1533963e+01 -4.5152342e+01 2.0027892e-01 0.0000000e+00 1.0000000e+01 2.8605457e+03 5.0486657e+00 0.0000000e+00 2.0027673e-01 1.1224508e+01 2.0647986e-01 3.3736235e-01 3.5979544e+02 1.6664060e-02 2.6785368e+02 -2.9934996e+01 5.0758087e+00 -3.1022000e+00 -2.5300458e-03 7.9249568e-03 0.0000000e+00 0 6 2 5 3.0 diff --git a/smoke_test/plotting.py b/smoke_test/plotting.py index 4aea79a..a05cf78 100644 --- a/smoke_test/plotting.py +++ b/smoke_test/plotting.py @@ -73,6 +73,78 @@ def _parse_header(lc_file: Path) -> Tuple[List[float] | None, List[float] | None return planet_vals, event_vals +def _sanity_check_flux_conservation( + lc_file: Path, + true_flux: np.ndarray, + src1_flux: np.ndarray | None, + src2_flux: np.ndarray | None, + abs_tol: float = 1e-5, +) -> None: + """Check that true_relative_flux ~= src1 + src2 + fblend for a couple of epochs. + + This is a strict, non-configurable sanity check. It only runs when both + source1/source2 columns exist and the header contains #fs and #fs2. On + failure it raises SmokeTestError so CI/smoke-test runner fails. + """ + # Only run when both source columns are present + if src1_flux is None or src2_flux is None: + return + + fs = None + fs2 = None + with lc_file.open(encoding="utf-8") as header_reader: + for raw in header_reader: + if not raw.startswith("#"): + break + s = raw.strip() + if s.startswith("#fs:"): + parts = s.split() + if len(parts) >= 2: + try: + fs = float(parts[1]) + except Exception: + fs = None + elif s.startswith("#fs2:"): + parts = s.split() + if len(parts) >= 2: + try: + fs2 = float(parts[1]) + except Exception: + fs2 = None + + if fs is None or fs2 is None: + # If header doesn't provide both, skip the strict check + return + + fblend = 1.0 - float(fs) - float(fs2) + + # choose representative epochs: first epoch and the epoch of peak true flux + try: + peak_idx = int(np.nanargmax(true_flux)) + except Exception: + peak_idx = 0 + + indices = [0] + if peak_idx != 0: + indices.append(int(peak_idx)) + + for idx in indices: + try: + expected = float(src1_flux[idx]) + float(src2_flux[idx]) + float(fblend) + truev = float(true_flux[idx]) + except Exception: + # If indexing or conversion fails, raise an error to surface the issue + raise SmokeTestError( + f"Flux sanity check failed: could not index/convert epoch {idx} in {lc_file.name}" + ) + diff = abs(truev - expected) + if diff > abs_tol: + raise SmokeTestError( + f"Flux sanity check failed for {lc_file.name} at epoch index {idx}: true={truev:.8f}, expected={expected:.8f}, abs_diff={diff:.8e}, abs_tol={abs_tol}" + ) + + + def _plot_photometry_only( lc_file: Path, output_dir: Path, @@ -81,6 +153,8 @@ def _plot_photometry_only( flux: np.ndarray, flux_err: np.ndarray, true_flux: np.ndarray | None, + src1_flux: np.ndarray | None = None, + src2_flux: np.ndarray | None = None, ) -> Path: fig, ax = plt.subplots(1, 1, figsize=(10, 6)) fig.suptitle(title, fontsize=14) @@ -106,6 +180,28 @@ def _plot_photometry_only( zorder=2, alpha=0.8, ) + if src1_flux is not None: + ax.plot( + time, + src1_flux, + "-", + linewidth=1.2, + color="magenta", + label="Source 1", + zorder=4, + alpha=0.7, + ) + if src2_flux is not None: + ax.plot( + time, + src2_flux, + "-", + linewidth=1.2, + color="gold", + label="Source 2", + zorder=5, + alpha=0.7, + ) ax.axhline(1.0, color="k", linestyle="--", linewidth=1.5, label="Baseline", zorder=3) ax.set_xlabel("Time (days)") ax.set_ylabel("Relative Flux") @@ -408,6 +504,8 @@ def _render_astrometric_figure( true_y_vals: np.ndarray | None, meas_x: np.ndarray | None, meas_y: np.ndarray | None, + src1_flux: np.ndarray | None = None, + src2_flux: np.ndarray | None = None, ) -> Tuple[Path, Path | None]: fig, axes = plt.subplots(2, 2, figsize=(12, 10)) fig.suptitle(title, fontsize=14) @@ -438,6 +536,28 @@ def _render_astrometric_figure( zorder=2, alpha=0.8, ) + if src1_flux is not None: + ax_light.plot( + time, + src1_flux, + "-", + linewidth=1.2, + color="magenta", + label="Source 1", + zorder=4, + alpha=0.7, + ) + if src2_flux is not None: + ax_light.plot( + time, + src2_flux, + "-", + linewidth=1.2, + color="gold", + label="Source 2", + zorder=5, + alpha=0.7, + ) ax_light.axhline( 1.0, color="k", @@ -695,6 +815,13 @@ def plot_lightcurves( if val is not None: astrometry_expected = str(val).strip().lower() not in {"0", "false", "off"} + # Determine whether MULTIPLE_SOURCES was requested in the parameter file. + multiple_sources_expected = False + if params: + ms = params.get("MULTIPLE_SOURCES") + if ms is not None: + multiple_sources_expected = str(ms).strip().lower() not in {"0", "false", "off"} + for lc_file in lc_files: planet_vals, event_vals = _parse_header(lc_file) df = pd.read_csv(lc_file, sep=r"\s+", comment="#") @@ -799,8 +926,23 @@ def _optional_column(name: str) -> np.ndarray | None: f"Smoke test failed: astrometric columns missing in {lc_file.name}" ) + # Extract source flux columns (optional for binary source events) + src1_flux = _optional_column("source1_relative_flux") + src2_flux = _optional_column("source2_relative_flux") + + # If the parameter file did not request multiple sources, ignore any + # source1/source2 columns that may nevertheless appear in the output + # (prevents plotting Source 1/Source 2 when MULTIPLE_SOURCES isn't set). + if not multiple_sources_expected: + src1_flux = None + src2_flux = None + + # Strict flux-conservation sanity check (non-configurable). Only run + # when multiple sources are expected and both columns are present. + _sanity_check_flux_conservation(lc_file, true_flux, src1_flux, src2_flux) + if not has_astrom: - _plot_photometry_only(lc_file, output_dir, title, time, flux, flux_err, true_flux) + _plot_photometry_only(lc_file, output_dir, title, time, flux, flux_err, true_flux, src1_flux, src2_flux) continue true_N_mas = _require_column("true_N_centroid_mas") @@ -945,6 +1087,8 @@ def _optional_column(name: str) -> np.ndarray | None: true_y_vals, meas_x, meas_y, + src1_flux, + src2_flux, ) if astrometry_expected: diff --git a/smoke_test/runner.py b/smoke_test/runner.py index 4e76f52..4631ae1 100644 --- a/smoke_test/runner.py +++ b/smoke_test/runner.py @@ -90,7 +90,9 @@ def _resolve_case_selection(raw_choices: Sequence[str] | None, ci_mode: bool = F ("smoke_std", "gulls_std.x", "smoke_std.prm"), ("smoke_std_binary", "gulls_std.x", "smoke_std_binary.prm"), ("smoke_fish", "gullsFish.x", "smoke_fish.prm"), - ("smoke_fish_binary", "gullsFish.x", "smoke_std_binary.prm"), + # For the fish binary CI case we should use the fish binary + # parameter file so outputs land under the fish/ output tree. + ("smoke_fish_binary", "gullsFish.x", "smoke_fish_binary.prm"), ("smoke_croin", "gulls_croin.x", "smoke_croin.prm"), ("smoke_croin_binary", "gulls_croin.x", "smoke_croin_binary.prm"), ) diff --git a/src/outputLightcurve.cpp b/src/outputLightcurve.cpp index 226079b..3444ee1 100644 --- a/src/outputLightcurve.cpp +++ b/src/outputLightcurve.cpp @@ -295,15 +295,14 @@ void outputLightcurve(struct event *Event, struct obsfilekeywords World[], struc { static const char* baseCols[] = { "Simulation_time", "measured_relative_flux", "measured_relative_flux_error", - "true_relative_flux", "true_relative_flux_error", "observatory_code", + "true_relative_flux", "true_relative_flux_error", "source1_relative_flux", "source2_relative_flux","observatory_code", "saturation_flag", "best_single_lens_fit", "x_centroid", "x_centroid_error","y_centroid", "y_centroid_error", "true_x_centroid", "true_x_centroid_error","true_y_centroid", "true_y_centroid_error", - "parallax_shift_t", - "parallax_shift_u", "BJD", "source_x", - "source_y", "source2_x", "source2_y", "lens1_x", "lens1_y", - "lens2_x", "lens2_y", "parallax_shift_x", - "parallax_shift_y", "parallax_shift_z" + "parallax_shift_t","parallax_shift_u", "BJD", + "source_x", "source_y", "source2_x", "source2_y", + "lens1_x", "lens1_y", "lens2_x", "lens2_y", + "parallax_shift_x", "parallax_shift_y", "parallax_shift_z" }; int nBase = sizeof(baseCols) / sizeof(baseCols[0]); for(int i = 0; i < nBase; ++i) { @@ -346,40 +345,68 @@ void outputLightcurve(struct event *Event, struct obsfilekeywords World[], struc } // finish the header line fprintf(lcfile_ptr, "\n"); - } - //output the lightcurve - int shiftedidx; + } + //output the lightcurve + int shiftedidx; - if(lcfile_ptr!=NULL && fileOpen==1) + if(lcfile_ptr!=NULL && fileOpen==1) + { + for(i=0;inepochs;i++) { - for(i=0;inepochs;i++) - { t=Event->epoch[i]; obsidx=Event->obsidx[i]; shiftedidx = i-Event->nepochsvec[obsidx]; - - fprintf(lcfile_ptr, "%.12g %.8g %g %.12g %g %d %d %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.6g %.6g %16.7f %.6g %.6g %.6g %.6g %.6g %.6g %.6g %.6g %.6g %.6g %.6g ", + + // compute per-source relative fluxes (fractions of baseline) + double src1_rel = 0.0, src2_rel = 0.0; + int sc_local = -1; + if(Event->scompanions.size()>0) sc_local = Event->scompanions[0]; + if(Paramfile->multiple_sources && sc_local>-1) + { + // filter for this epoch + int filt = World[obsidx].filter; + double r = 0.0; // because it's zero is there's no companion? + // if the event structure has the flux ratio and it's the right size for the number of filters + if(Event->scomp_fsofs1.size()>0 && Event->scomp_fsofs1[0].size()>filt) + // r is the flux ratio for the current filter + r = Event->scomp_fsofs1[0][filt]; //fs2/fs1 + double FS1 = Event->fs[obsidx]; // is this fs1 or fs1+fs2? It's fs1 + double FS2 = FS1 * r; + // Use stored magnifications from lightcurve generation + src1_rel = FS1 * Event->Asrc1[i]; + src2_rel = FS2 * Event->Asrc2[i]; + } + fprintf(lcfile_ptr, + "%.12g %.8g %g " + "%.12g %g %.8g %.8g %d" + "%d %.8g " + "%.8g %.8g %.8g %.8g " + "%.8g %.8g %.8g %.8g " + "%.6g %.6g %16.7f " + "%.6g %.6g %.6g %.6g " + "%.6g %.6g %.6g %.6g " + "%.6g %.6g %.6g ", Event->epoch[i], Event->Aobs[i], Event->Aerr[i], //0, 1, 2 - Event->Atrue[i], Event->Atrueerr[i], obsidx, //3, 4, 5 - (Event->nosat[i]?0:1), Event->Afit[i], //6, 7 - Event->xc[i], Event->xcerr[i], Event->yc[i], Event->ycerr[i], - Event->xctrue[i], Event->xctrueerr[i], Event->yctrue[i], Event->yctrueerr[i], + Event->Atrue[i], Event->Atrueerr[i], src1_rel, src2_rel, obsidx, //3, 4, 5, 6, 7 + (Event->nosat[i]?0:1), Event->Afit[i], //8, 9 + Event->xc[i], Event->xcerr[i], Event->yc[i], Event->ycerr[i], //10, 11, 12, 13 + Event->xctrue[i], Event->xctrueerr[i], Event->yctrue[i], Event->yctrueerr[i], //14, 15, 16, 17 //Event->pllx[obsidx].tshift(Event->jdepoch[i]), //8 //Event->pllx[obsidx].ushift(Event->jdepoch[i]), //9 //Event->pllx[obsidx].epochs[Event->jdepoch[i]], //10 - Event->pllx[obsidx].tshift[shiftedidx], - Event->pllx[obsidx].ushift[shiftedidx], - Event->pllx[obsidx].epochs[shiftedidx], - Event->xs[i], Event->ys[i], - Event->xs2[i], Event->ys2[i], - Event->xl1[i], Event->yl1[i], - Event->xl2[i], Event->yl2[i], //14, 15, 16 + Event->pllx[obsidx].tshift[shiftedidx], //18 + Event->pllx[obsidx].ushift[shiftedidx], //19 + Event->pllx[obsidx].epochs[shiftedidx], //20 + Event->xs[i], Event->ys[i], //21, 22 + Event->xs2[i], Event->ys2[i], //23, 24 + Event->xl1[i], Event->yl1[i], //25, 26 + Event->xl2[i], Event->yl2[i], //27, 28 //Event->pllx[obsidx].sslocation[Event->jdepoch[i]][0], //17 //Event->pllx[obsidx].sslocation[Event->jdepoch[i]][1], //18 //Event->pllx[obsidx].sslocation[Event->jdepoch[i]][2]); //19 - Event->pllx[obsidx].sslocation[shiftedidx][0], //17 - Event->pllx[obsidx].sslocation[shiftedidx][1], //18 - Event->pllx[obsidx].sslocation[shiftedidx][2]); //19 + Event->pllx[obsidx].sslocation[shiftedidx][0], //29 + Event->pllx[obsidx].sslocation[shiftedidx][1], //30 + Event->pllx[obsidx].sslocation[shiftedidx][2]); //31 if(ndF>0) @@ -464,9 +491,11 @@ void outputImages(struct event *Event, struct obsfilekeywords World[], struct sl } else { - tmp1 = "%s%s_%d_%d_%d",Paramfile->outputdir + Paramfile->run_name + "_" - + to_string(Event->instance) + "_" + to_string(Paramfile->choosefield) + "_" - + to_string(Event->id); + // Build the same base name as the other branch — avoid stray format literal and + // do not use the comma operator. Keep it as a plain concatenation. + tmp1 = Paramfile->outputdir + Paramfile->run_name + "_" + + to_string(Event->instance) + "_" + to_string(Paramfile->choosefield) + "_" + + to_string(Event->id); } basefname=tmp1; @@ -474,13 +503,14 @@ void outputImages(struct event *Event, struct obsfilekeywords World[], struct sl for(int obsidx=0;obsidxnumobservatories;obsidx++) { if(Event->nepochsvec[obsidx+1]-Event->nepochsvec[obsidx]<=0) continue; - tmp1 += "." + to_string(obsidx) + "_"; + // do not mutate tmp1/basefname; create a small suffix for this observation + string obs_suffix = "." + to_string(obsidx) + "_"; filter = World[obsidx].filter; //first the baseline image imtype="base"; - oname = basefname + tmp1 + imtype + extension + ".fits"; + oname = basefname + obs_suffix + imtype + extension + ".fits"; mag = Sources->mags[Event->source][filter]; diff --git a/src/pllxLightcurveGenerator.cpp b/src/pllxLightcurveGenerator.cpp index 884e904..98885b2 100644 --- a/src/pllxLightcurveGenerator.cpp +++ b/src/pllxLightcurveGenerator.cpp @@ -41,9 +41,9 @@ void lightcurveGenerator(struct filekeywords* Paramfile, struct event *Event, st //if the event is saturated in each band, no need to calculate the lightcurve if(Event->nepochs==0 || Event->allsat) - { - return; - } + { + return; + } vector idxshift; int shiftedidx; const double timeout = Paramfile->lc_timeout; @@ -55,16 +55,16 @@ void lightcurveGenerator(struct filekeywords* Paramfile, struct event *Event, st //Calculate the lightcurve for (int idx = 0; idx < Event->nepochs; ++idx) - { - if (enforce_timeout) - { - time_t now = time(NULL); - if (difftime(now, starttime) > timeout) - { - timed_out = true; - break; - } - } + { + if (enforce_timeout) + { + time_t now = time(NULL); + if (difftime(now, starttime) > timeout) + { + timed_out = true; + break; + } + } obsidx = Event->obsidx[idx]; shiftedidx=idx-idxshift[obsidx]; double amp = 0.0; @@ -72,28 +72,27 @@ void lightcurveGenerator(struct filekeywords* Paramfile, struct event *Event, st { // Lightcurve is identical from observatory to observatory amp = Event->Atrue[shiftedidx]; - } - else - { + } else { + double tt = (Event->epoch[idx] - Event->t0) / Event->tE_r; double uu = u0; if (Paramfile->pllxMultiplyer) - { + { tt += Event->pllx[obsidx].tshift[shiftedidx]; uu += Event->pllx[obsidx].ushift[shiftedidx]; - } + } Event->umin=min(Event->umin,qAdd(tt,uu)); double xsCoM = tt * cosa - uu * sina + VBM_origin; - double ysCenter = tt * sina + uu * cosa; + double ysCenter = tt * sina + uu * cosa; - Event->xs[idx] = xsCoM; - Event->ys[idx] = ysCenter; - Event->xl1[idx] = VBM_origin; - Event->yl1[idx] = 0.0; - Event->xl2[idx] = VBM_origin + a; - Event->yl2[idx] = 0.0; + Event->xs[idx] = xsCoM; + Event->ys[idx] = ysCenter; + Event->xl1[idx] = VBM_origin; + Event->yl1[idx] = 0.0; + Event->xl2[idx] = VBM_origin + a; + Event->yl2[idx] = 0.0; Event->vbm->a1 = lim_gamma; amp = Event->vbm->BinaryMag2(a, q, xsCoM, ysCenter, rs); @@ -102,42 +101,47 @@ void lightcurveGenerator(struct filekeywords* Paramfile, struct event *Event, st Event->vbm_therr[idx] = Event->vbm->therr; if(Paramfile->multiple_sources && Event->scompanions.size()>0) - { - double xs2CoM, ys2Center; - double x2off = Event->scomp_s[0] * cos(Event->scomp_phase[0]*TO_RAD); - double y2off = Event->scomp_s[0] * sin(Event->scomp_phase[0]*TO_RAD) * cos(Event->scomp_inc[0]*TO_RAD); - xs2CoM = xsCoM + x2off * cos(Event->scomp_alpha[0]*TO_RAD) - y2off * sin(Event->scomp_alpha[0]*TO_RAD); - ys2Center = ysCenter + x2off * sin(Event->scomp_alpha[0]*TO_RAD) + y2off * cos(Event->scomp_alpha[0]*TO_RAD); + { + double xs2CoM, ys2Center; + double x2off = Event->scomp_s[0] * cos(Event->scomp_phase[0]*TO_RAD); + double y2off = Event->scomp_s[0] * sin(Event->scomp_phase[0]*TO_RAD) * cos(Event->scomp_inc[0]*TO_RAD); + xs2CoM = xsCoM + x2off * cos(Event->scomp_alpha[0]*TO_RAD) - y2off * sin(Event->scomp_alpha[0]*TO_RAD); + ys2Center = ysCenter + x2off * sin(Event->scomp_alpha[0]*TO_RAD) + y2off * cos(Event->scomp_alpha[0]*TO_RAD); - double amp2 = Event->vbm->BinaryMag2(a, q, xs2CoM, ys2Center, Event->scomp_rs[0]); - Event->xs2[idx] = xs2CoM; - Event->ys2[idx] = ys2Center; - - int filt = World[obsidx].filter; - Event->Atrue[idx] = amp + Event->scomp_fsofs1[0][filt] * (amp2-1); + double amp2 = Event->vbm->BinaryMag2(a, q, xs2CoM, ys2Center, Event->scomp_rs[0]); + Event->xs2[idx] = xs2CoM; + Event->ys2[idx] = ys2Center; - //Put binary source astrometry here - //Event->xctrue[idx] = ; - //Event->yctrue[idx] = ; + // Store individual source magnifications + Event->Asrc1[idx] = amp; + Event->Asrc2[idx] = amp2; + + int filt = World[obsidx].filter; + Event->Atrue[idx] = amp + Event->scomp_fsofs1[0][filt] * (amp2-1); - } - else - { - Event->Atrue[idx] = amp; + //Put binary source astrometry here + //Event->xctrue[idx] = ; + //Event->yctrue[idx] = ; + } else { + + Event->Asrc1[idx] = amp; + Event->Asrc2[idx] = 0.0; // No second source + Event->Atrue[idx] = amp; //put single source astrometry here //Event->xctrue[idx] = ; //Event->yctrue[idx] = ; - } - } + + } + } - // Keep track of highest magnification - if (amp > Event->Amax) { - Event->Amax = amp; - Event->peakpoint = idx; - } + // Keep track of highest magnification + if (amp > Event->Amax) { + Event->Amax = amp; + Event->peakpoint = idx; } +} if(timed_out) { diff --git a/src/structures.h b/src/structures.h index 77fc97d..b29d9a9 100644 --- a/src/structures.h +++ b/src/structures.h @@ -358,6 +358,8 @@ struct event{ vector Aobs; vector Aerr; vector Afit; + vector Asrc1; //magnification of source 1 + vector Asrc2; //magnification of source 2 vector nosat; /*Is point unsaturated? */ vector backmag; vector dF; diff --git a/src/timeSequencer.cpp b/src/timeSequencer.cpp index b53ef45..115010a 100644 --- a/src/timeSequencer.cpp +++ b/src/timeSequencer.cpp @@ -259,6 +259,8 @@ void setupMemory(struct obsfilekeywords World[], struct event *Event, struct fil Event->yl1.resize(Event->nepochs); Event->xl2.resize(Event->nepochs); Event->yl2.resize(Event->nepochs); + Event->Asrc1.resize(Event->nepochs); + Event->Asrc2.resize(Event->nepochs); Event->vbm_rootaccuracy.resize(Event->nepochs); Event->vbm_squarecheck.resize(Event->nepochs); Event->vbm_therr.resize(Event->nepochs);