From 44db3d4d648d58f70c25abc57a1ea202e43fc32b Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 10 Sep 2025 18:43:51 +0200 Subject: [PATCH 01/40] Add files via upload Elementi aggiunti in row 16/17 double tanomaly; double resanomaly; --- RTModel/include/bumper.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RTModel/include/bumper.h b/RTModel/include/bumper.h index 8fd70e2..bc02778 100644 --- a/RTModel/include/bumper.h +++ b/RTModel/include/bumper.h @@ -13,6 +13,8 @@ class bumper{ double *dp; double *curv, *cov; double Amp; + double tanomaly; + double resanomaly; int nps; char modelcode[16]; char *buffer; From 6e33cfa6f77b6c6dfaf62c73056329dc0f74156a Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 10 Sep 2025 18:45:46 +0200 Subject: [PATCH 02/40] Add files via upload Elementi aggiunti in row 46/47 double tmaxmax; double maxmaxsum; --- RTModel/include/LevMarFit.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RTModel/include/LevMarFit.h b/RTModel/include/LevMarFit.h index 8f4125d..1b3fa2c 100644 --- a/RTModel/include/LevMarFit.h +++ b/RTModel/include/LevMarFit.h @@ -43,6 +43,8 @@ class LevMar { double* constraints, * consleft, * consright, * consvars; int modnumber; + double tmaxmax; + double maxmaxsum; double Tol; // void (LevMar::* PrintOut)(double*); From 4d17a277d13ffc2953b9179e3a128a8d694f46bd Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 10 Sep 2025 18:47:09 +0200 Subject: [PATCH 03/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 83 ++++++++++++++++++++++++++++++++------- 1 file changed, 69 insertions(+), 14 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 9237509..9ae497d 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -11,7 +11,6 @@ #include #include #include -#include #include #include @@ -239,11 +238,12 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[3] = log(pr[3]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); } + else { modnumber = 0; nps = 4; @@ -258,11 +258,12 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[3] = log(pr[3]); } current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); break; + case 'B': if (modelcode[1] == 'O') { modnumber = 3; @@ -277,8 +278,9 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[0] = log(pr[0]); pr[1] = log(pr[1]); pr[6] = log(pr[6]); + current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -298,7 +300,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[6] = log(pr[6]); } current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); @@ -338,7 +340,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[4] = log(pr[4]); pr[5] = log(pr[5]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -358,7 +360,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[4] = log(pr[4]); pr[5] = log(pr[5]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -397,7 +399,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[7] = log(pr[7]); pr[8] = log(pr[8]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -1303,6 +1305,10 @@ void LevMar::EvaluateModel(double* pr, int fl, int ips) { double LevMar::ChiSquared(double* pr) { double chi2 = 0, chi0, chia, p1; + double p1max = 0, maxsum = 0; + double tmax = 0; + + maxmaxsum = 0; for (int fl = 0; fl < nfil; fl++) { VBM->satellite = satel[starts[fl]]; @@ -1376,12 +1382,57 @@ double LevMar::ChiSquared(double* pr) { p1 = (y[i] - pr[nps + filter[i] * nlinpar] - pr[nps + 1 + filter[i] * nlinpar] * fb[i]) * w[i] * w[i] * pr[nps + 1 + filter[i] * nlinpar] * Tol; chi0 += p1 * p1; p1 = (y[i] - pr[nps + filter[i] * nlinpar] - pr[nps + 1 + filter[i] * nlinpar] * fb[i]) * w[i]; - chi2 += p1 * p1; - if (pr[nps + 1 + filter[i] * nlinpar] > 2 * y[i]) { - flagblending++; - chi2 += (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * w[i] * w[i]; + + if (p1 > 0) { + maxsum += p1; // somma i residui + + if (p1 > p1max) { + p1max = p1; // massimo residuo positivo + //std::cout << "\nIl massimo residuo positivo è: " << p1max << " al punto " << i << std::endl; + tmax = t[i]; // tempo in cui si ha il massimo residuo positivo + } + + } + else { + //printf("\nResiduo negativo o nullo %d: %lf", i, p1); // stampa il residuo negativo + // Reset in caso di interruzione della sequenza positiva + maxsum = 0; // azzera la somma dei residui + } + //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente + if (maxsum > maxmaxsum) { + maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi + //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; + tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + } + //std::cout << std::endl; } + //else { + //printf("\nNessun residuo positivo trovato."); + //} + + + //robustezza del fit + //bisogna considerare il residuo di una seguenza di punti consecutivi che hanno lo stesso segno + //**** il segno è dato dal confronto tra il modello ed il flusso misurato + //**** se il modello è minore del flusso misurato allora il flusso ha un picco che il modello non ha + //**** viceversa si ha un deep + //quindi bisogna considerare la somma dei residui dei punti consecutivi che hanno lo stesso segno + //in questo modo bisogna tener conto del segno di p1 ad ogni passo e confrontarlo con il segno del residuo precedente + //se ha lo stesso segno allora si sommano i residui + // in caso di più sequenze si considera qualla in cui la somma dei residui è max + //infine si calcola la massima deviazione standard dei residui all'interno di una sequenza di punti consecutivi che hanno lo stesso segno + + chi2 += p1 * p1; + if (pr[nps + 1 + filter[i] * nlinpar] > 2 * y[i]) { + flagblending++; + chi2 += (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * w[i] * w[i]; + } + } + if (maxsum > maxmaxsum) { + maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi + //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; + tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } chi0 = sqrt(2 * chi0); // Error in chi square if (chi0 / chi2 > 0.1) Tol *= 0.5; @@ -1752,6 +1803,11 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { pr[i] = ((pr[i] > -1.e300) && (pr[i] < 1.e300)) ? pr[i] : -1.e300; fprintf(f, "%le ", pr[i]); } + //stampa tmaxmax + fprintf(f, "%.16le ", tmaxmax); + + //stampa maxmaxsum + fprintf(f, "%.16le ", maxmaxsum); // Write chi square fprintf(f, "%.16le\n", c0); @@ -1792,5 +1848,4 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { fprintf(f, "\n"); } fclose(f); -} - +} \ No newline at end of file From 47b3c47d73260be863f92386d27aec51eaa8f8f2 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 10 Sep 2025 18:50:49 +0200 Subject: [PATCH 04/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index 847721f..f4fe47a 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -7,7 +7,6 @@ #include #include #include -#include #include #include #include "bumper.h" @@ -26,7 +25,7 @@ int main(int argc, char* argv[]) { char modelcode[256]; char command[256], buffer[10000]; double value; - double t, y, w, errDecv, * pr, * sigmapr, * Cov, * Curv, fac, facr, c1, c0, chithr, bestplan = 1.e100, bestbin = 1.e100; + double t, y, w, errDecv, * pr, * sigmapr, * Cov, * Curv, fac, facr, c1, tmaxmax, maxmaxsum, c0, chithr, bestplan = 1.e100, bestbin = 1.e100; double renorm; int nfil, il, nlc, nmod, np, k; int nps = 4, dof, nlinpar; @@ -201,6 +200,10 @@ int main(int argc, char* argv[]) { for (int j = 0; j < nps + nlinpar * nfil; j++) { fscanf(f, "%le", &(pr[j])); } + + //per leggere i parametri + fscanf(f, "%le", &(tmaxmax)); + fscanf(f, "%le", &(maxmaxsum)); fscanf(f, "%le", &(c0)); renorm = sqrt(c0 / dof); @@ -273,6 +276,8 @@ int main(int argc, char* argv[]) { bumperlist->SetCovariance(Cov, dof / c0); strcpy(bumperlist->modelcode, filename); bumperlist->SetBuffer(f, start, end); + bumperlist->tanomaly = tmaxmax; + bumperlist->maxsum = maxmaxsum; bumperlist->Amp = c0; bumperlist->next = scanbumper; } @@ -283,6 +288,8 @@ int main(int argc, char* argv[]) { scanbumper2->SetCovariance(Cov, dof / c0); strcpy(scanbumper2->modelcode, filename); scanbumper2->SetBuffer(f, start, end); + bumperlist->tanomaly = tmaxmax; + bumperlist->maxsum = maxmaxsum; scanbumper2->Amp = c0; scanbumper2->next = scanbumper->next; scanbumper->next = scanbumper2; @@ -293,6 +300,8 @@ int main(int argc, char* argv[]) { bumperlist->SetCovariance(Cov, dof / c0); strcpy(bumperlist->modelcode, filename); bumperlist->SetBuffer(f, start, end); + bumperlist->tanomaly = tmaxmax; + bumperlist->maxsum = maxmaxsum; bumperlist->Amp = c0; } nmod++; From 673d1221944f2c9efb5d3833460c1c09c3e3280b Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 10 Sep 2025 18:53:10 +0200 Subject: [PATCH 05/40] Update LevMarFit.cpp --- RTModel/lib/LevMarFit.cpp | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 9ae497d..61be3e6 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -238,7 +238,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[3] = log(pr[3]); current_path(exedir); - //current_path(".."); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -258,7 +258,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[3] = log(pr[3]); } current_path(exedir); - //current_path(".."); + current_path(".."); current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); @@ -280,7 +280,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[6] = log(pr[6]); current_path(exedir); - //current_path(".."); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -300,7 +300,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[6] = log(pr[6]); } current_path(exedir); - //current_path(".."); + current_path(".."); current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); @@ -340,7 +340,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[4] = log(pr[4]); pr[5] = log(pr[5]); current_path(exedir); - //current_path(".."); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -360,7 +360,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[4] = log(pr[4]); pr[5] = log(pr[5]); current_path(exedir); - //current_path(".."); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -399,7 +399,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[7] = log(pr[7]); pr[8] = log(pr[8]); current_path(exedir); - //current_path(".."); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -1388,7 +1388,7 @@ double LevMar::ChiSquared(double* pr) { if (p1 > p1max) { p1max = p1; // massimo residuo positivo - //std::cout << "\nIl massimo residuo positivo è: " << p1max << " al punto " << i << std::endl; + //std::cout << "\nIl massimo residuo positivo è: " << p1max << " al punto " << i << std::endl; tmax = t[i]; // tempo in cui si ha il massimo residuo positivo } @@ -1399,10 +1399,10 @@ double LevMar::ChiSquared(double* pr) { maxsum = 0; // azzera la somma dei residui } - //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente + //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente if (maxsum > maxmaxsum) { maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi - //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; + //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } //std::cout << std::endl; @@ -1414,13 +1414,13 @@ double LevMar::ChiSquared(double* pr) { //robustezza del fit //bisogna considerare il residuo di una seguenza di punti consecutivi che hanno lo stesso segno - //**** il segno è dato dal confronto tra il modello ed il flusso misurato - //**** se il modello è minore del flusso misurato allora il flusso ha un picco che il modello non ha + //**** il segno è dato dal confronto tra il modello ed il flusso misurato + //**** se il modello è minore del flusso misurato allora il flusso ha un picco che il modello non ha //**** viceversa si ha un deep //quindi bisogna considerare la somma dei residui dei punti consecutivi che hanno lo stesso segno //in questo modo bisogna tener conto del segno di p1 ad ogni passo e confrontarlo con il segno del residuo precedente //se ha lo stesso segno allora si sommano i residui - // in caso di più sequenze si considera qualla in cui la somma dei residui è max + // in caso di più sequenze si considera qualla in cui la somma dei residui è max //infine si calcola la massima deviazione standard dei residui all'interno di una sequenza di punti consecutivi che hanno lo stesso segno chi2 += p1 * p1; @@ -1431,7 +1431,7 @@ double LevMar::ChiSquared(double* pr) { } if (maxsum > maxmaxsum) { maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi - //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; + //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } chi0 = sqrt(2 * chi0); // Error in chi square @@ -1848,4 +1848,5 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { fprintf(f, "\n"); } fclose(f); -} \ No newline at end of file + +} From 71de4493f9627fefc8abccd4c73cf1ee87df10fe Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 10 Sep 2025 19:00:11 +0200 Subject: [PATCH 06/40] Update LevMarFit.cpp --- RTModel/lib/LevMarFit.cpp | 38 +++++++++----------------------------- 1 file changed, 9 insertions(+), 29 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 61be3e6..7ed7c25 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1388,51 +1388,30 @@ double LevMar::ChiSquared(double* pr) { if (p1 > p1max) { p1max = p1; // massimo residuo positivo - //std::cout << "\nIl massimo residuo positivo è: " << p1max << " al punto " << i << std::endl; tmax = t[i]; // tempo in cui si ha il massimo residuo positivo } } else { - //printf("\nResiduo negativo o nullo %d: %lf", i, p1); // stampa il residuo negativo - // Reset in caso di interruzione della sequenza positiva - maxsum = 0; // azzera la somma dei residui + maxsum = 0; } - //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente if (maxsum > maxmaxsum) { maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi - //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - //std::cout << std::endl; + } - //else { - //printf("\nNessun residuo positivo trovato."); - //} - - - //robustezza del fit - //bisogna considerare il residuo di una seguenza di punti consecutivi che hanno lo stesso segno - //**** il segno è dato dal confronto tra il modello ed il flusso misurato - //**** se il modello è minore del flusso misurato allora il flusso ha un picco che il modello non ha - //**** viceversa si ha un deep - //quindi bisogna considerare la somma dei residui dei punti consecutivi che hanno lo stesso segno - //in questo modo bisogna tener conto del segno di p1 ad ogni passo e confrontarlo con il segno del residuo precedente - //se ha lo stesso segno allora si sommano i residui - // in caso di più sequenze si considera qualla in cui la somma dei residui è max - //infine si calcola la massima deviazione standard dei residui all'interno di una sequenza di punti consecutivi che hanno lo stesso segno - - chi2 += p1 * p1; + + chi2 += p1 * p1; if (pr[nps + 1 + filter[i] * nlinpar] > 2 * y[i]) { flagblending++; chi2 += (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * w[i] * w[i]; } } if (maxsum > maxmaxsum) { - maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi - //std::cout << "\nLa somma massima dei residui positivi consecutivi è: " << maxmaxsum << " al tempo " << tmax << std::endl; - tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + maxmaxsum = maxsum; + tmaxmax = tmax; } chi0 = sqrt(2 * chi0); // Error in chi square if (chi0 / chi2 > 0.1) Tol *= 0.5; @@ -1803,10 +1782,10 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { pr[i] = ((pr[i] > -1.e300) && (pr[i] < 1.e300)) ? pr[i] : -1.e300; fprintf(f, "%le ", pr[i]); } - //stampa tmaxmax + // Write tmaxmax fprintf(f, "%.16le ", tmaxmax); - //stampa maxmaxsum + // Write maxmaxsum fprintf(f, "%.16le ", maxmaxsum); // Write chi square fprintf(f, "%.16le\n", c0); @@ -1850,3 +1829,4 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { fclose(f); } + From 024a1dae08e341040735af81dea15f489fccf9c2 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Tue, 16 Sep 2025 10:08:52 +0200 Subject: [PATCH 07/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 36 +++++++++++++++++++----------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 7ed7c25..eb22fc7 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -238,7 +238,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[3] = log(pr[3]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -258,7 +258,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[3] = log(pr[3]); } current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); @@ -280,7 +280,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[6] = log(pr[6]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -300,7 +300,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[6] = log(pr[6]); } current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); @@ -340,7 +340,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[4] = log(pr[4]); pr[5] = log(pr[5]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -360,7 +360,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[4] = log(pr[4]); pr[5] = log(pr[5]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -399,7 +399,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[7] = log(pr[7]); pr[8] = log(pr[8]); current_path(exedir); - current_path(".."); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -1301,6 +1301,7 @@ void LevMar::EvaluateModel(double* pr, int fl, int ips) { } break; } + } double LevMar::ChiSquared(double* pr) { @@ -1393,25 +1394,28 @@ double LevMar::ChiSquared(double* pr) { } else { - maxsum = 0; + maxsum = 0; // azzera la somma dei residui } + //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente if (maxsum > maxmaxsum) { - maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi + maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } } + + - chi2 += p1 * p1; + chi2 += p1 * p1; if (pr[nps + 1 + filter[i] * nlinpar] > 2 * y[i]) { flagblending++; chi2 += (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * w[i] * w[i]; } } if (maxsum > maxmaxsum) { - maxmaxsum = maxsum; - tmaxmax = tmax; + maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi + tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } chi0 = sqrt(2 * chi0); // Error in chi square if (chi0 / chi2 > 0.1) Tol *= 0.5; @@ -1782,10 +1786,10 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { pr[i] = ((pr[i] > -1.e300) && (pr[i] < 1.e300)) ? pr[i] : -1.e300; fprintf(f, "%le ", pr[i]); } - // Write tmaxmax + //stampa tmaxmax fprintf(f, "%.16le ", tmaxmax); - // Write maxmaxsum + //stampa maxmaxsum fprintf(f, "%.16le ", maxmaxsum); // Write chi square fprintf(f, "%.16le\n", c0); @@ -1827,6 +1831,4 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { fprintf(f, "\n"); } fclose(f); - -} - +} \ No newline at end of file From fcb7bf41719546d38c83db6e14e281d96b223620 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 24 Sep 2025 12:05:28 +0200 Subject: [PATCH 08/40] Add files via upload From 4fe212a5379d00570343a10d636fc1fbaa0b800c Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 24 Sep 2025 12:14:27 +0200 Subject: [PATCH 09/40] Add files via upload --- RTModel/include/bumper.h | 1 + 1 file changed, 1 insertion(+) diff --git a/RTModel/include/bumper.h b/RTModel/include/bumper.h index bc02778..295f5f5 100644 --- a/RTModel/include/bumper.h +++ b/RTModel/include/bumper.h @@ -15,6 +15,7 @@ class bumper{ double Amp; double tanomaly; double resanomaly; + double maxsum, maxmaxsum, tmaxmax, tmax; int nps; char modelcode[16]; char *buffer; From 94bfdc954ab1cc45d9e8df8a3d622b3cbfb3dffe Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Tue, 4 Nov 2025 17:24:22 +0100 Subject: [PATCH 10/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 65 ++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 31 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index f4fe47a..ed10335 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -43,7 +43,7 @@ int main(int argc, char* argv[]) { // Directory preliminaries. Reads event name from arguments. - auto exedir = current_path(); + auto exedir = std::filesystem::current_path(); if (argc > 2) { strcpy(eventname, argv[1]); @@ -129,7 +129,6 @@ int main(int argc, char* argv[]) { } } fclose(f); - // Determine model type nps = (astrometric) ? 4 : 0; nlinpar = (astrometric) ? 4 : 2; @@ -165,6 +164,7 @@ int main(int argc, char* argv[]) { supfac *= nps; + printf("\n- Model code: %s", modelcode); @@ -808,22 +808,23 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } - double mindpeak = 1.e100; - int idpeak = 0; - for (int ipeak = 0; ipeak < npeaks; ipeak++) { - double ddpeak = fabs(peaks[ipeak] - pr[2]); - if (ddpeak < mindpeak) { - idpeak = ipeak; - mindpeak = ddpeak; - } - } + //double mindpeak = 1.e100; + //int idpeak = 0; + //for (int ipeak = 0; ipeak < npeaks; ipeak++) { + // double ddpeak = fabs(peaks[ipeak] - pr[2]); + // if (ddpeak < mindpeak) { + // idpeak = ipeak; + // mindpeak = ddpeak; + // } + //} double u0 = exp(pr[0]); double s0, s; - for (int ipeak = 0; ipeak < npeaks; ipeak++) { - if (ipeak == idpeak) continue; + /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { + if (ipeak == idpeak) continue;*/ double q = 0.001; - double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); + double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); + /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ double xc0 = sqrt(u0 * u0 + dt * dt), xc; double alpha0 = atan2(u0, -dt), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); @@ -856,7 +857,7 @@ int main(int argc, char* argv[]) { fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[1]), pr[2]); alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[1]), pr[2]); - } + //} scanbumper = scanbumper->next; } fclose(g); @@ -892,7 +893,7 @@ int main(int argc, char* argv[]) { } fclose(f); - double mindpeak = 1.e100; + /*double mindpeak = 1.e100; int idpeak = 0; for (int ipeak = 0; ipeak < npeaks; ipeak++) { double ddpeak = fabs(peaks[ipeak] - pr[2]); @@ -900,7 +901,7 @@ int main(int argc, char* argv[]) { idpeak = ipeak; mindpeak = ddpeak; } - } + }*/ scanbumper = bumperlist; for (il = 1; il <= nmod; il++) { @@ -908,7 +909,7 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } - double mindpeak = 1.e100; + /*double mindpeak = 1.e100; int idpeak = 0; for (int ipeak = 0; ipeak < npeaks; ipeak++) { double ddpeak = fabs(peaks[ipeak] - pr[2]); @@ -916,14 +917,15 @@ int main(int argc, char* argv[]) { idpeak = ipeak; mindpeak = ddpeak; } - } + }*/ double u0 = pr[0]; double s0, s; - for (int ipeak = 0; ipeak < npeaks; ipeak++) { - if (ipeak == idpeak) continue; + /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { + if (ipeak == idpeak) continue;*/ double q = 0.001; - double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); + double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); + //double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); double xc0 = sqrt(u0 * u0 + dt * dt), xc; double alpha0 = atan2(u0, -dt), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); @@ -972,7 +974,7 @@ int main(int argc, char* argv[]) { if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); fprintf(g, "\n"); - } + //} scanbumper = scanbumper->next; } @@ -1007,7 +1009,7 @@ int main(int argc, char* argv[]) { } fclose(f); - double mindpeak = 1.e100; + /*double mindpeak = 1.e100; int idpeak = 0; for (int ipeak = 0; ipeak < npeaks; ipeak++) { double ddpeak = fabs(peaks[ipeak] - pr[2]); @@ -1015,7 +1017,7 @@ int main(int argc, char* argv[]) { idpeak = ipeak; mindpeak = ddpeak; } - } + }*/ scanbumper = bumperlist; for (il = 1; il <= nmod; il++) { @@ -1023,7 +1025,7 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } - double mindpeak = 1.e100; + /*double mindpeak = 1.e100; int idpeak = 0; for (int ipeak = 0; ipeak < npeaks; ipeak++) { double ddpeak = fabs(peaks[ipeak] - pr[2]); @@ -1031,14 +1033,15 @@ int main(int argc, char* argv[]) { idpeak = ipeak; mindpeak = ddpeak; } - } + }*/ double u0 = pr[0]; double s0, s; - for (int ipeak = 0; ipeak < npeaks; ipeak++) { - if (ipeak == idpeak) continue; + /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { + if (ipeak == idpeak) continue;*/ double q = 0.001; - double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); + double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); + //double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); double xc0 = sqrt(u0 * u0 + dt * dt), xc; double alpha0 = atan2(u0, -dt), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); @@ -1088,7 +1091,7 @@ int main(int argc, char* argv[]) { if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); fprintf(g, "\n"); - } + //} scanbumper = scanbumper->next; } From 49f5692c624fcc0911a99483f0c3bf3efb155b92 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 12 Nov 2025 13:55:19 +0100 Subject: [PATCH 11/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 170 +++++++++++++++++++++++++--------- 1 file changed, 128 insertions(+), 42 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index ed10335..1fe9710 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -14,6 +14,72 @@ using namespace std; using namespace std::filesystem; +static inline void initcond_from_anomaly( + FILE* g, + double xc0, + double q, + double u0, + double alpha0, + double rho, + double* pr, + int nps, + int astrometric +) { + double xc = xc0; + double s0 = 0.5 * (sqrt(4 + xc * xc) + xc); + while (xc < 4 * sqrt(q) / (s0 * s0)) q *= 0.1; + + //prima parte + xc = xc0 + 4 * sqrt(q) / (s0 * s0); + double s = 0.5 * (sqrt(4 + xc * xc) + xc); + double alpha = alpha0; + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + if (astrometric) + fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + fprintf(g, "\n"); + + //seconda parte + xc = xc0 + 4 * sqrt(q) / (s0 * s0); + s = 0.5 * (sqrt(4 + xc * xc) + xc); + alpha = alpha0; + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + if (astrometric) + fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + fprintf(g, "\n"); + + xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) - xc); + q = 0.001; + while (xc < 3 * sqrt(3 * q) * s0 * s0 * s0) q *= 0.1; + + //terza parte + xc = xc0 - 3 * sqrt(3 * q) * s0 * s0 * s0; + s = 0.5 * (sqrt(4 + xc * xc) - xc); + alpha = alpha0 + M_PI + asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + fprintf(g, "\n"); + + alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + fprintf(g, "\n"); + + //quarta parte + xc = xc0 + 3 * sqrt(3 * q) * s0 * s0 * s0; + s = 0.5 * (sqrt(4 + xc * xc) - xc); + alpha = alpha0 + M_PI + asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + fprintf(g, "\n"); + + alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + fprintf(g, "\n"); +} + double supfac = 3.0; // factor multiplying the inverse covariance in search for superpositions (models are incompatible if farther than supfac*sigma) double chifac = 1.; // number of sigmas in the chi square distribution for accepting alternative models after the best one int maxmodels = 10; // maximum number of models returned @@ -129,6 +195,7 @@ int main(int argc, char* argv[]) { } } fclose(f); + // Determine model type nps = (astrometric) ? 4 : 0; nlinpar = (astrometric) ? 4 : 2; @@ -168,7 +235,6 @@ int main(int argc, char* argv[]) { printf("\n- Model code: %s", modelcode); - pr = (double*)malloc(sizeof(double) * (nps + 20 + nlinpar * nfil)); //20 added to allow for higher order models in updates sigmapr = (double*)malloc(sizeof(double) * (nps + nlinpar * nfil)); Cov = (double*)malloc(sizeof(double) * (nps * nps)); @@ -185,6 +251,7 @@ int main(int argc, char* argv[]) { return 0; } current_path("PreModels"); + nmod = 0; auto searchstring = regex(string(modelcode) + ".*txt"); @@ -314,6 +381,7 @@ int main(int argc, char* argv[]) { printf("\n! No models for this class"); return 0; } + printf("\nnmod: %d", nmod); c1 = (bumperlist) ? bumperlist->Amp : 1.e100; printf("\nNumber of initial models = %d\nMinimum chi square = %lf\n", nmod, c1); @@ -521,7 +589,7 @@ int main(int argc, char* argv[]) { if (nmod > maxmodels) nmod = maxmodels; printf("\nModels to be saved = %d\n", nmod); - + // Store best models passing the selection in directory "Models" @@ -704,6 +772,7 @@ int main(int argc, char* argv[]) { } } } + if (modelcode[1] == 'O') { printf("\n- Preparing initial conditions for Keplerian orbital motion"); if (f = fopen("InitCondLK.txt", "r")) { @@ -930,49 +999,63 @@ int main(int argc, char* argv[]) { double alpha0 = atan2(u0, -dt), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); - xc = xc0; - s0 = 0.5 * (sqrt(4 + xc * xc) + xc); - while (xc < 4 * sqrt(q) / (s0 * s0)) q *= 0.1; - xc = xc0 + 4 * sqrt(q) / (s0 * s0); - s = 0.5 * (sqrt(4 + xc * xc) + xc); - alpha = alpha0; - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); - if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); - fprintf(g, "\n"); - xc = xc0 - 4 * sqrt(q) / (s0 * s0); - s = 0.5 * (sqrt(4 + xc * xc) + xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); - if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); - fprintf(g, "\n"); + //funzione inline + initcond_from_anomaly(g, xc0, q, u0, alpha0, rho, pr, nps, astrometric); + //xc = xc0; - xc = xc0; - s0 = 0.5 * (sqrt(4 + xc * xc) - xc); - q = 0.001; - while (xc < 3 * sqrt(3 * q) * s0 * s0 * s0) q *= 0.1; - xc = xc0 - 3 * sqrt(3 * q) * s0 * s0 * s0; - s = 0.5 * (sqrt(4 + xc * xc) - xc); - alpha = alpha0 + M_PI + asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); - if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); - fprintf(g, "\n"); - alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); - if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); - fprintf(g, "\n"); + //s0 = 0.5 * (sqrt(4 + xc * xc) + xc); + //while (xc < 4 * sqrt(q) / (s0 * s0)) q *= 0.1; - xc = xc0 + 3 * sqrt(3 * q) * s0 * s0 * s0; - s = 0.5 * (sqrt(4 + xc * xc) - xc); - alpha = alpha0 + M_PI + asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); - if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); - fprintf(g, "\n"); + //xc = xc0 + 4 * sqrt(q) / (s0 * s0); + //s = 0.5 * (sqrt(4 + xc * xc) + xc); + //alpha = alpha0; + //fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + //if (astrometric) + // fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + //fprintf(g, "\n"); - alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); - if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); - fprintf(g, "\n"); + ////seconda parte + //xc = xc0 - 4 * sqrt(q) / (s0 * s0); + //s = 0.5 * (sqrt(4 + xc * xc) + xc); + //fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + //if (astrometric) + // fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + //fprintf(g, "\n"); + + + //xc = xc0; + + //s0 = 0.5 * (sqrt(4 + xc * xc) - xc); + //q = 0.001; + //while (xc < 3 * sqrt(3 * q) * s0 * s0 * s0) q *= 0.1; + + ////terza parte + //xc = xc0 - 3 * sqrt(3 * q) * s0 * s0 * s0; + //s = 0.5 * (sqrt(4 + xc * xc) - xc); + //alpha = alpha0 + M_PI + asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + //fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + //if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + //fprintf(g, "\n"); + + //alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + //fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + //if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + //fprintf(g, "\n"); + + ////quarta parte + //xc = xc0 + 3 * sqrt(3 * q) * s0 * s0 * s0; + //s = 0.5 * (sqrt(4 + xc * xc) - xc); + //alpha = alpha0 + M_PI + asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + //fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + //if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + //fprintf(g, "\n"); + + //alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); + //fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); + //if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + //fprintf(g, "\n"); //} scanbumper = scanbumper->next; @@ -1046,7 +1129,10 @@ int main(int argc, char* argv[]) { double alpha0 = atan2(u0, -dt), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); - xc = xc0; + //funzione inline + initcond_from_anomaly(g, xc0, q, u0, alpha0, rho, pr, nps, astrometric); + + /* xc = xc0; s0 = 0.5 * (sqrt(4 + xc * xc) + xc); while (xc < 4 * sqrt(q) / (s0 * s0)) q *= 0.1; @@ -1089,7 +1175,7 @@ int main(int argc, char* argv[]) { alpha = alpha0 + M_PI - asin(fabs(2 * sqrt(q * (1 - s * s)) / s) / xc); fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le 0.0001 0.0001 0.0001", s, q, u0, alpha, rho, exp(pr[1]), pr[2], pr[4], pr[5]); if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); - fprintf(g, "\n"); + fprintf(g, "\n");*/ //} scanbumper = scanbumper->next; From f2b08c13eaf4f1d25b0dc803858fd8209c0723c8 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Sun, 16 Nov 2025 12:13:05 +0100 Subject: [PATCH 12/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 83 ++++++++++++++++++++++++++++----------- 1 file changed, 59 insertions(+), 24 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index eb22fc7..fba7af4 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -5,7 +5,7 @@ #define _USE_MATH_DEFINES #include "LevMarFit.h" #include "bumper.h" -#include +#include "VBMicrolensingLibrary.h" #include #include #include @@ -188,6 +188,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { printf("\n\n- Event: %s\n", eventname); current_path(eventname); + printf("- Working directory: %s\n", current_path().string().c_str()); //if(argc>2){ // strcpy(eventname,argv[1]); @@ -202,7 +203,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { current_path(eventname); current_path("Data"); - + printf("- Data directory: %s\n", current_path().string().c_str()); /* Reading coordinates */ auto searchstring = regex(".*\\.coordinates"); @@ -214,8 +215,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { break; } } - current_path(eventname); + printf("\n- Event directory: %s\n", current_path().string().c_str()); // Read light curve ReadCurve(); @@ -237,11 +238,12 @@ void LevMar::ReadFiles(int argc, char* argv[]) { error = InitCond(presigmapr, preleftlim, prerightlim); pr[1] = log(pr[1]); pr[3] = log(pr[3]); - current_path(exedir); + //current_path(exedir); //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); + } else { @@ -257,11 +259,13 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[3] = log(pr[3]); } - current_path(exedir); + //current_path(exedir); //current_path(".."); + current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); + break; case 'B': @@ -279,7 +283,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[6] = log(pr[6]); - current_path(exedir); + //current_path(exedir); //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); @@ -299,7 +303,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[6] = log(pr[6]); } - current_path(exedir); + //current_path(exedir); //current_path(".."); current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); @@ -319,8 +323,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[4] = log(pr[4]); pr[5] = log(pr[5]); - current_path(exedir); - current_path(".."); + //current_path(exedir); + //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -339,7 +343,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[4] = log(pr[4]); pr[5] = log(pr[5]); - current_path(exedir); + //current_path(exedir); //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); @@ -359,7 +363,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[4] = log(pr[4]); pr[5] = log(pr[5]); - current_path(exedir); + //current_path(exedir); //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); @@ -398,11 +402,12 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[5] = log(pr[5]); pr[7] = log(pr[7]); pr[8] = log(pr[8]); - current_path(exedir); + //current_path(exedir); //current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); + } else { modnumber = 8; @@ -1289,6 +1294,9 @@ void LevMar::EvaluateModel(double* pr, int fl, int ips) { VBM->BinaryLightCurveKepler(pr, tfl, fbfl, y1a, y2a, seps, sizes[fl]); } break; + + + case 8: VBM->TripleLightCurve(pr, tfl, fbfl, y1a, y2a, sizes[fl]); break; @@ -1306,7 +1314,8 @@ void LevMar::EvaluateModel(double* pr, int fl, int ips) { double LevMar::ChiSquared(double* pr) { double chi2 = 0, chi0, chia, p1; - double p1max = 0, maxsum = 0; + double p1maxp = 0, maxsump = 0; + double p1maxn = 0, maxsumn = 0; double tmax = 0; maxmaxsum = 0; @@ -1384,25 +1393,42 @@ double LevMar::ChiSquared(double* pr) { chi0 += p1 * p1; p1 = (y[i] - pr[nps + filter[i] * nlinpar] - pr[nps + 1 + filter[i] * nlinpar] * fb[i]) * w[i]; + if (p1 > 0) { - maxsum += p1; // somma i residui + maxsumn = 0; // azzera la somma dei residui negativi + maxsump += p1; // somma i residui positivi - if (p1 > p1max) { - p1max = p1; // massimo residuo positivo + if (p1 > p1maxp) { + p1maxp = p1; // aggiorna il massimo residuo positivo tmax = t[i]; // tempo in cui si ha il massimo residuo positivo } } else { - maxsum = 0; // azzera la somma dei residui + maxsump = 0; // azzera la somma dei residui positivi + maxsumn += p1; // somma dei residui negativi + if (abs(p1) > abs(p1maxn)) { + p1maxn = p1; // aggiorna il massimo residuo negativo + + //salva il tempo del primo e del'ultimo elemento della sequenza + + } } //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente - if (maxsum > maxmaxsum) { - maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi - tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + if (maxsump > maxmaxsum) { + if (maxsump > abs(maxsumn)) { + maxmaxsum = maxsump; + tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + } + } + else (abs(maxsumn) > maxmaxsum); { + if (abs(maxsumn) > maxsump) { + maxmaxsum = maxsumn; + //calcolo del tempo + + } } - } @@ -1413,9 +1439,18 @@ double LevMar::ChiSquared(double* pr) { chi2 += (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * (pr[nps + 1 + filter[i] * nlinpar] - 2 * y[i]) * w[i] * w[i]; } } - if (maxsum > maxmaxsum) { - maxmaxsum = maxsum; // somma massima dei residui positivi consecutivi - tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + if (maxsump > maxmaxsum) { + if (maxsump > abs(maxsumn)) { + maxmaxsum = maxsump; + tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + } + } + else (abs(maxsumn) > maxmaxsum); { + if (abs(maxsumn) > maxsump) { + maxmaxsum = maxsumn; + //calcolo del tempo + + } } chi0 = sqrt(2 * chi0); // Error in chi square if (chi0 / chi2 > 0.1) Tol *= 0.5; From 134dc0d9b7eebdc644f0f8188c2e3de97c6b5f96 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Tue, 16 Dec 2025 09:23:59 +0100 Subject: [PATCH 13/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 83 +++++++++++++++++++++++++-------------- 1 file changed, 54 insertions(+), 29 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index fba7af4..c568d7c 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1314,9 +1314,9 @@ void LevMar::EvaluateModel(double* pr, int fl, int ips) { double LevMar::ChiSquared(double* pr) { double chi2 = 0, chi0, chia, p1; - double p1maxp = 0, maxsump = 0; - double p1maxn = 0, maxsumn = 0; - double tmax = 0; + double p1max = 0, maxsump = 0; + double maxsumn = 0; + double t1 = 0, t2 = 0, tmax = 0; maxmaxsum = 0; @@ -1386,6 +1386,8 @@ double LevMar::ChiSquared(double* pr) { } // Photometric chi square + bool in_pos_sequence; + flagblending = 0; for (int i = 0; i < np; i++) { if (w[i] > 0) { @@ -1393,46 +1395,66 @@ double LevMar::ChiSquared(double* pr) { chi0 += p1 * p1; p1 = (y[i] - pr[nps + filter[i] * nlinpar] - pr[nps + 1 + filter[i] * nlinpar] * fb[i]) * w[i]; + //controlli su i e filter[i] + //azzerare tutti gli elementi maxsumn e maxsump + //inizializzare t1 a t[i] + //inizializzare in_pos_sequence = (p1>0) + if ((i == 0) || (filter[i] != filter[i-1])) { + maxsumn = 0; + maxsump = 0; + t1 = t[i]; + } if (p1 > 0) { + if (!in_pos_sequence) { + in_pos_sequence = true; + + t2 = t[i]; // tempo primo positivo dopo una sequenza neg + + } maxsumn = 0; // azzera la somma dei residui negativi maxsump += p1; // somma i residui positivi - if (p1 > p1maxp) { - p1maxp = p1; // aggiorna il massimo residuo positivo + if (p1 > p1max) { //per trovare il picco + p1max = p1; // aggiorna il massimo residuo positivo tmax = t[i]; // tempo in cui si ha il massimo residuo positivo } } else { - maxsump = 0; // azzera la somma dei residui positivi - maxsumn += p1; // somma dei residui negativi + //fine sequenza pos + if (in_pos_sequence) { + in_pos_sequence = false; + t1 = t[i-1]; // tempo ultimo pos prima di una sequenza neg + } + //aggiorna t2 + t2 = t[i]; - if (abs(p1) > abs(p1maxn)) { - p1maxn = p1; // aggiorna il massimo residuo negativo + maxsump = 0; // azzera la somma dei residui positivi + maxsumn += p1; // somma dei residui negativi, è un numero negativo + p1max = 0; - //salva il tempo del primo e del'ultimo elemento della sequenza - - } } //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente if (maxsump > maxmaxsum) { - if (maxsump > abs(maxsumn)) { - maxmaxsum = maxsump; - tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi - } + maxmaxsum = maxsump; + tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + } else (abs(maxsumn) > maxmaxsum); { - if (abs(maxsumn) > maxsump) { - maxmaxsum = maxsumn; - //calcolo del tempo - + maxmaxsum = - maxsumn; + //calcolo del tempo || pr[6] è t0 in binary lens + if (abs(pr[6] - t1) > abs(pr[6] - t2)) { + tmaxmax = t1; + } + else { + tmaxmax = t2; } } } + pr[2]; - chi2 += p1 * p1; if (pr[nps + 1 + filter[i] * nlinpar] > 2 * y[i]) { flagblending++; @@ -1440,18 +1462,20 @@ double LevMar::ChiSquared(double* pr) { } } if (maxsump > maxmaxsum) { - if (maxsump > abs(maxsumn)) { - maxmaxsum = maxsump; - tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi - } + maxmaxsum = maxsump; + tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi + } else (abs(maxsumn) > maxmaxsum); { - if (abs(maxsumn) > maxsump) { - maxmaxsum = maxsumn; - //calcolo del tempo - + maxmaxsum = maxsumn; + if (abs(pr[6] - t1) > abs(pr[6] - t2)) { + tmaxmax = t1; + } + else { + tmaxmax = t2; } } + chi0 = sqrt(2 * chi0); // Error in chi square if (chi0 / chi2 > 0.1) Tol *= 0.5; if (chi0 / chi2 < 0.01 && Tol < .99e-2) Tol *= 2; @@ -1826,6 +1850,7 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { //stampa maxmaxsum fprintf(f, "%.16le ", maxmaxsum); + // Write chi square fprintf(f, "%.16le\n", c0); From da8892d69caf8584e219ce13556f143243c8b243 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Tue, 16 Dec 2025 10:08:41 +0100 Subject: [PATCH 14/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index c568d7c..e9c3ab3 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1441,10 +1441,10 @@ double LevMar::ChiSquared(double* pr) { tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - else (abs(maxsumn) > maxmaxsum); { + else if(fabs(maxsumn) > maxmaxsum); { maxmaxsum = - maxsumn; //calcolo del tempo || pr[6] è t0 in binary lens - if (abs(pr[6] - t1) > abs(pr[6] - t2)) { + if (abs(double(pr[6] - t1)) > abs(double(pr[6] - t2))) { tmaxmax = t1; } else { @@ -1466,7 +1466,7 @@ double LevMar::ChiSquared(double* pr) { tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - else (abs(maxsumn) > maxmaxsum); { + else (fabs(double(maxsumn)) > maxmaxsum); { maxmaxsum = maxsumn; if (abs(pr[6] - t1) > abs(pr[6] - t2)) { tmaxmax = t1; From 2aec2712029d84167cf698d1766d74b145735a8f Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Thu, 18 Dec 2025 12:41:48 +0100 Subject: [PATCH 15/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index e9c3ab3..0df2965 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1466,8 +1466,8 @@ double LevMar::ChiSquared(double* pr) { tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - else (fabs(double(maxsumn)) > maxmaxsum); { - maxmaxsum = maxsumn; + else if (fabs(maxsumn) > maxmaxsum); { + maxmaxsum = -maxsumn; if (abs(pr[6] - t1) > abs(pr[6] - t2)) { tmaxmax = t1; } From 04c6d159f850ed6d5e6003848207840cea45d7e3 Mon Sep 17 00:00:00 2001 From: Valerio Bozza Date: Thu, 18 Dec 2025 14:54:23 +0100 Subject: [PATCH 16/40] Update pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 5ae48e6..e31afdd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,7 +11,7 @@ authors = [ { name = "Valerio Bozza", email = "valboz@sa.infn.it" }, ] license = { text = "GPL-3.0" } -requires-python = ">=3.7,<4" +requires-python = ">=3.8,<4" readme = "README.md" classifiers = [ 'Development Status :: 5 - Production/Stable', From 2793e4456cbd6b0e63eeb395b53ac28304412d65 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Sun, 21 Dec 2025 15:12:29 +0100 Subject: [PATCH 17/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 0df2965..100b3ed 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1453,7 +1453,6 @@ double LevMar::ChiSquared(double* pr) { } } - pr[2]; chi2 += p1 * p1; if (pr[nps + 1 + filter[i] * nlinpar] > 2 * y[i]) { From 130906e466b7da0507ab5d22810a74e23cf87dcc Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Sun, 21 Dec 2025 17:02:14 +0100 Subject: [PATCH 18/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 66 +++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 100b3ed..fb51e6f 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1406,13 +1406,13 @@ double LevMar::ChiSquared(double* pr) { } if (p1 > 0) { - if (!in_pos_sequence) { - in_pos_sequence = true; + //if (!in_pos_sequence) { + // in_pos_sequence = true; - t2 = t[i]; // tempo primo positivo dopo una sequenza neg - - } - maxsumn = 0; // azzera la somma dei residui negativi + // t2 = t[i]; // tempo primo positivo dopo una sequenza neg + // + //} + /*maxsumn = 0; */ // azzera la somma dei residui negativi maxsump += p1; // somma i residui positivi if (p1 > p1max) { //per trovare il picco @@ -1423,16 +1423,16 @@ double LevMar::ChiSquared(double* pr) { } else { //fine sequenza pos - if (in_pos_sequence) { - in_pos_sequence = false; - t1 = t[i-1]; // tempo ultimo pos prima di una sequenza neg - } - //aggiorna t2 - t2 = t[i]; + //if (in_pos_sequence) { + // in_pos_sequence = false; + // t1 = t[i-1]; // tempo ultimo pos prima di una sequenza neg + //} + ////aggiorna t2 + //t2 = t[i]; maxsump = 0; // azzera la somma dei residui positivi - maxsumn += p1; // somma dei residui negativi, è un numero negativo - p1max = 0; + //maxsumn += p1; // somma dei residui negativi, è un numero negativo + //p1max = 0; } //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente @@ -1441,16 +1441,16 @@ double LevMar::ChiSquared(double* pr) { tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - else if(fabs(maxsumn) > maxmaxsum); { - maxmaxsum = - maxsumn; - //calcolo del tempo || pr[6] è t0 in binary lens - if (abs(double(pr[6] - t1)) > abs(double(pr[6] - t2))) { - tmaxmax = t1; - } - else { - tmaxmax = t2; - } - } + //else if(fabs(maxsumn) > maxmaxsum); { + // maxmaxsum = - maxsumn; + // //calcolo del tempo || pr[6] è t0 in binary lens + // if (abs(double(pr[6] - t1)) > abs(double(pr[6] - t2))) { + // tmaxmax = t1; + // } + // else { + // tmaxmax = t2; + // } + //} } @@ -1465,15 +1465,15 @@ double LevMar::ChiSquared(double* pr) { tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - else if (fabs(maxsumn) > maxmaxsum); { - maxmaxsum = -maxsumn; - if (abs(pr[6] - t1) > abs(pr[6] - t2)) { - tmaxmax = t1; - } - else { - tmaxmax = t2; - } - } + //else if (fabs(maxsumn) > maxmaxsum); { + // maxmaxsum = -maxsumn; + // if (abs(pr[6] - t1) > abs(pr[6] - t2)) { + // tmaxmax = t1; + // } + // else { + // tmaxmax = t2; + // } + //} chi0 = sqrt(2 * chi0); // Error in chi square if (chi0 / chi2 > 0.1) Tol *= 0.5; From 930bbb1a776004bd5eabd10f8b442ebe26a63aa1 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Sun, 21 Dec 2025 22:58:29 +0100 Subject: [PATCH 19/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index fb51e6f..48c69b8 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1399,11 +1399,11 @@ double LevMar::ChiSquared(double* pr) { //azzerare tutti gli elementi maxsumn e maxsump //inizializzare t1 a t[i] //inizializzare in_pos_sequence = (p1>0) - if ((i == 0) || (filter[i] != filter[i-1])) { + /*if ((i == 0) || (filter[i] != filter[i-1])) { maxsumn = 0; maxsump = 0; t1 = t[i]; - } + }*/ if (p1 > 0) { //if (!in_pos_sequence) { @@ -1466,7 +1466,7 @@ double LevMar::ChiSquared(double* pr) { } //else if (fabs(maxsumn) > maxmaxsum); { - // maxmaxsum = -maxsumn; + // maxmaxsum = - maxsumn; // if (abs(pr[6] - t1) > abs(pr[6] - t2)) { // tmaxmax = t1; // } From daf50337f6886c06f59d892c14b97afa61ac91b2 Mon Sep 17 00:00:00 2001 From: Valerio Bozza Date: Fri, 26 Dec 2025 17:30:56 +0100 Subject: [PATCH 20/40] Update build_wheels.yml --- .github/workflows/build_wheels.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build_wheels.yml b/.github/workflows/build_wheels.yml index 6cb9a8c..0b157ed 100644 --- a/.github/workflows/build_wheels.yml +++ b/.github/workflows/build_wheels.yml @@ -15,7 +15,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [ubuntu-latest, windows-latest, macos-13,macos-14] # was macos-latest + os: [ubuntu-latest, windows-latest, macos-14,macos-15] # was macos-latest steps: - uses: actions/checkout@v4 From 4405fc75606c3978799f8db2ff4c624518c8a806 Mon Sep 17 00:00:00 2001 From: Valerio Bozza Date: Fri, 26 Dec 2025 17:48:23 +0100 Subject: [PATCH 21/40] Update build_wheels.yml --- .github/workflows/build_wheels.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build_wheels.yml b/.github/workflows/build_wheels.yml index 0b157ed..e4cf42e 100644 --- a/.github/workflows/build_wheels.yml +++ b/.github/workflows/build_wheels.yml @@ -15,7 +15,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [ubuntu-latest, windows-latest, macos-14,macos-15] # was macos-latest + os: [ubuntu-latest, windows-latest,macos-latest] # was macos-latest steps: - uses: actions/checkout@v4 From 50689faf02b0864851ec0b6fc362bb85979acaff Mon Sep 17 00:00:00 2001 From: Valerio Bozza Date: Fri, 26 Dec 2025 17:58:03 +0100 Subject: [PATCH 22/40] Update pyproject.toml --- pyproject.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index e31afdd..1c77c37 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,5 +38,8 @@ minimum-version = "build-system.requires" [tool.cibuildwheel] build-verbosity = "3" +[tool.cibuildwheel.macos] +archs = ["x86_64", "arm64"] + [tool.setuptools.package-data] "RTModel.data" = ["data/*"] From e600e3b5cf384bed2ffe4265c96add8784224a30 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Fri, 26 Dec 2025 19:26:24 +0100 Subject: [PATCH 23/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 68 +++++++++++++++++++-------------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 48c69b8..1595b83 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1399,20 +1399,20 @@ double LevMar::ChiSquared(double* pr) { //azzerare tutti gli elementi maxsumn e maxsump //inizializzare t1 a t[i] //inizializzare in_pos_sequence = (p1>0) - /*if ((i == 0) || (filter[i] != filter[i-1])) { + if ((i == 0) || (filter[i] != filter[i-1])) { maxsumn = 0; maxsump = 0; t1 = t[i]; - }*/ + } if (p1 > 0) { - //if (!in_pos_sequence) { - // in_pos_sequence = true; + if (!in_pos_sequence) { + in_pos_sequence = true; - // t2 = t[i]; // tempo primo positivo dopo una sequenza neg - // - //} - /*maxsumn = 0; */ // azzera la somma dei residui negativi + t2 = t[i]; // tempo primo positivo dopo una sequenza neg + + } + maxsumn = 0; // azzera la somma dei residui negativi maxsump += p1; // somma i residui positivi if (p1 > p1max) { //per trovare il picco @@ -1423,16 +1423,16 @@ double LevMar::ChiSquared(double* pr) { } else { //fine sequenza pos - //if (in_pos_sequence) { - // in_pos_sequence = false; - // t1 = t[i-1]; // tempo ultimo pos prima di una sequenza neg - //} + if (in_pos_sequence) { + in_pos_sequence = false; + t1 = t[i-1]; // tempo ultimo pos prima di una sequenza neg + } ////aggiorna t2 - //t2 = t[i]; + t2 = t[i]; maxsump = 0; // azzera la somma dei residui positivi - //maxsumn += p1; // somma dei residui negativi, è un numero negativo - //p1max = 0; + maxsumn += p1; // somma dei residui negativi, è un numero negativo + p1max = 0; } //Alla fine del ciclo, se la somma dei residui positivi consecutivi è maggiore della somma massima trovata finora allora aggiorna la somma massima e il tempo corrispondente @@ -1441,16 +1441,16 @@ double LevMar::ChiSquared(double* pr) { tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - //else if(fabs(maxsumn) > maxmaxsum); { - // maxmaxsum = - maxsumn; - // //calcolo del tempo || pr[6] è t0 in binary lens - // if (abs(double(pr[6] - t1)) > abs(double(pr[6] - t2))) { - // tmaxmax = t1; - // } - // else { - // tmaxmax = t2; - // } - //} + else if(fabs(maxsumn) > maxmaxsum); { + maxmaxsum = - maxsumn; + //calcolo del tempo || pr[6] è t0 in binary lens + if (abs(double(pr[6] - t1)) > abs(double(pr[6] - t2))) { + tmaxmax = t1; + } + else { + tmaxmax = t2; + } + } } @@ -1465,15 +1465,15 @@ double LevMar::ChiSquared(double* pr) { tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi } - //else if (fabs(maxsumn) > maxmaxsum); { - // maxmaxsum = - maxsumn; - // if (abs(pr[6] - t1) > abs(pr[6] - t2)) { - // tmaxmax = t1; - // } - // else { - // tmaxmax = t2; - // } - //} + else if (fabs(maxsumn) > maxmaxsum); { + maxmaxsum = - maxsumn; + if (abs(pr[6] - t1) > abs(pr[6] - t2)) { + tmaxmax = t1; + } + else { + tmaxmax = t2; + } + } chi0 = sqrt(2 * chi0); // Error in chi square if (chi0 / chi2 > 0.1) Tol *= 0.5; From f5f8c17cbb98206a0d2bc80a2d9a19bddbce4cf7 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Sun, 18 Jan 2026 11:52:01 +0100 Subject: [PATCH 24/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 255 ++++++++++++++++++++++++++++++---- 1 file changed, 231 insertions(+), 24 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index 1fe9710..9616b48 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -91,7 +91,8 @@ int main(int argc, char* argv[]) { char modelcode[256]; char command[256], buffer[10000]; double value; - double t, y, w, errDecv, * pr, * sigmapr, * Cov, * Curv, fac, facr, c1, tmaxmax, maxmaxsum, c0, chithr, bestplan = 1.e100, bestbin = 1.e100; + double t, y, w, errDecv, * pr, * sigmapr, * Cov, * Curv, fac, facr, c1, chithr, bestplan = 1.e100, bestbin = 1.e100; + double tmaxmax, y1maxmax, y2maxmax, maxmaxsum, c0; double renorm; int nfil, il, nlc, nmod, np, k; int nps = 4, dof, nlinpar; @@ -270,6 +271,8 @@ int main(int argc, char* argv[]) { //per leggere i parametri fscanf(f, "%le", &(tmaxmax)); + fscanf(f, "%le", &(y1maxmax)); + fscanf(f, "%le", &(y2maxmax)); fscanf(f, "%le", &(maxmaxsum)); fscanf(f, "%le", &(c0)); @@ -332,6 +335,12 @@ int main(int argc, char* argv[]) { //} } break; + + //da vedere per il triple lens + case 'T': + + + break; } @@ -344,6 +353,8 @@ int main(int argc, char* argv[]) { strcpy(bumperlist->modelcode, filename); bumperlist->SetBuffer(f, start, end); bumperlist->tanomaly = tmaxmax; + bumperlist->y1anomaly = y1maxmax; + bumperlist->y2anomaly = y2maxmax; bumperlist->maxsum = maxmaxsum; bumperlist->Amp = c0; bumperlist->next = scanbumper; @@ -356,6 +367,8 @@ int main(int argc, char* argv[]) { strcpy(scanbumper2->modelcode, filename); scanbumper2->SetBuffer(f, start, end); bumperlist->tanomaly = tmaxmax; + bumperlist->y1anomaly = y1maxmax; + bumperlist->y2anomaly = y2maxmax; bumperlist->maxsum = maxmaxsum; scanbumper2->Amp = c0; scanbumper2->next = scanbumper->next; @@ -368,6 +381,8 @@ int main(int argc, char* argv[]) { strcpy(bumperlist->modelcode, filename); bumperlist->SetBuffer(f, start, end); bumperlist->tanomaly = tmaxmax; + bumperlist->y1anomaly = y1maxmax; + bumperlist->y2anomaly = y2maxmax; bumperlist->maxsum = maxmaxsum; bumperlist->Amp = c0; } @@ -489,7 +504,7 @@ int main(int argc, char* argv[]) { if (facr < fac) fac = facr; scanbumper->signCovariance(1); break; - // Binary source + // Triple lens case 'T': // Check with no reflections pr[1] = scanbumper->p0[1]; @@ -507,8 +522,12 @@ int main(int argc, char* argv[]) { } } fac = facr; + // Need to implement all mass permutations (and reflections for static model) + break; + + //binary source case 'B': // Check with no reflections pr[1] = scanbumper->p0[1]; @@ -708,6 +727,106 @@ int main(int argc, char* argv[]) { } } } + printf("\n- Adding initial conditions for planets"); + if (f = fopen("InitCondTS.txt", "r")) { + int npeaks = 0; + double* peaks; + g = fopen("InitCondTS-temp.txt", "w"); + fscanf(f, "%d %d", &npeaks, &np); + if (npeaks == 0) { + fclose(g); + fclose(f); + remove("InitCondTS-temp.txt"); + } + else { + fprintf(g, "%d %d\n", npeaks, np + ((npeaks > 0) ? 6 * nmod : 0)); + peaks = (double*)malloc(sizeof(double) * npeaks); + + printf("\nNumber of initial conditions: %d", np + ((npeaks > 0) ? 6 * nmod : 0)); + for (int i = 0; i < npeaks; i++) { + fscanf(f, "%lg", &peaks[i]); + fprintf(g, "%le", peaks[i]); + fscanf(f, "%[^\n]s", &buffer); + fprintf(g, "%s\n", buffer); + } + for (int i = 0; i < np; i++) { + for (int j = 0; j < 10; j++) { + fscanf(f, "%lg", &pr[j]); + fprintf(g, "%.10le ", pr[j]); + } + fprintf(g, "\n"); + } + fclose(f); + + + scanbumper = bumperlist; + for (il = 1; il <= nmod; il++) { + for (int i = 0; i < nps; i++) { + pr[i] = scanbumper->p0[i]; + } + + double u0 = pr[2]; + double alpha = pr[3]; + double s0, s = exp(pr[0]), s2; + double q = exp(pr[1]); + double q2 = 0.000001; + /*double dt = (scanbumper->tanomaly - pr[6]) / exp(pr[5]);*/ + /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ + double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); + double xc0, xc; + double rho = exp(pr[4]); + double y1A, y2A; + double beta, beta0; + + //y1A e y2A posizione a tanomaly + //mentre y1anomaly e y2anomaly posizione ricavate dal LevMarfit + y1A = pr[2] * sin(pr[3]) - dt * cos(pr[3]) + q * s / (1 + q); //rispetto la primaria + y2A = -pr[2] * cos(pr[3]) - dt * sin(pr[3]); + beta = beta0 = atan2(y2A, y1A); + + xc0 = sqrt(y1A * y1A + y2A * y2A); + xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) + xc); + while (xc < 4 * sqrt(q2) / (s0 * s0)) q2 *= 0.1; + xc = xc0 + 4 * sqrt(q2) / (s0 * s0); + s2 = 0.5 * (sqrt(4 + xc * xc) + xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + xc = xc0 - 4 * sqrt(q2) / (s0 * s0); + s2 = 0.5 * (sqrt(4 + xc * xc) + xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + + xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) - xc); + q2 = 0.00001; + while (xc < 3 * sqrt(3 * q2) * s0 * s0 * s0) q2 *= 0.1; + xc = xc0 - 3 * sqrt(3 * q2) * s0 * s0 * s0; + s2 = 0.5 * (sqrt(4 + xc * xc) - xc); + beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + xc = xc0 + 3 * sqrt(3 * q2) * s0 * s0 * s0; + s2 = 0.5 * (sqrt(4 + xc * xc) - xc); + beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + + scanbumper = scanbumper->next; + } + fclose(g); + free(peaks); + remove("InitCondTS.txt"); + rename("InitCondTS-temp.txt", "InitCondTS.txt"); + } + } + if (modelcode[1] == 'X') { printf("\n- Preparing initial conditions for orbital motion"); if (f = fopen("InitCondLO.txt", "r")) { @@ -717,7 +836,7 @@ int main(int argc, char* argv[]) { fprintf(g, "%d %d\n", npeaks, np + nmod); printf("\nNumber of initial conditions: %d", np + nmod); for (int i = 0; i < np; i++) { - for (int j = 0; j < nps+3; j++) { + for (int j = 0; j < nps + 3; j++) { fscanf(f, "%lg", &pr[j]); fprintf(g, "%.10le ", pr[j]); } @@ -731,7 +850,7 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le 0.0e-006 0.0e-006 1.0e-006", exp(pr[0]), exp(pr[1]), pr[2], pr[3], exp(pr[4]), exp(pr[5]), pr[6], pr[7], pr[8]); - if(astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); fprintf(g, "\n"); scanbumper = scanbumper->next; } @@ -747,7 +866,7 @@ int main(int argc, char* argv[]) { fprintf(g, "%d %d\n", npeaks, np + nmod); printf("\nNumber of initial conditions: %d", np + nmod); for (int i = 0; i < np; i++) { - for (int j = 0; j < nps+5; j++) { + for (int j = 0; j < nps + 5; j++) { fscanf(f, "%lg", &pr[j]); fprintf(g, "%.10le ", pr[j]); } @@ -761,7 +880,7 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le 0.0e-006 0.0e-006 1.0e-006 1.0 1.0", exp(pr[0]), exp(pr[1]), pr[2], pr[3], exp(pr[4]), exp(pr[5]), pr[6], pr[7], pr[8]); - if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps-4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); + if (astrometric) fprintf(g, " %.10le %.10le %.10le %.10le", pr[nps - 4], pr[nps - 3], pr[nps - 2], pr[nps - 1]); fprintf(g, "\n"); scanbumper = scanbumper->next; @@ -771,8 +890,101 @@ int main(int argc, char* argv[]) { rename("InitCondLK-temp.txt", "InitCondLK.txt"); } } - } + ////bisogna aggiungere TX + + printf("\n- Adding initial conditions for planets"); + if (f = fopen("InitCondTX.txt", "r")) { + int npeaks = 0; + double* peaks; + g = fopen("InitCondTX-temp.txt", "w"); + fscanf(f, "%d %d", &npeaks, &np); + fprintf(g, "%d %d\n", npeaks, np + ((npeaks > 0) ? 6 * nmod : 0)); + peaks = (double*)malloc(sizeof(double) * npeaks); + + printf("\nNumber of initial conditions: %d", np + ((npeaks > 0) ? 6 * nmod : 0)); + for (int i = 0; i < npeaks; i++) { + fscanf(f, "%lg", &peaks[i]); + fprintf(g, "%le", peaks[i]); + fscanf(f, "%[^\n]s", &buffer); + fprintf(g, "%s\n", buffer); + } + for (int i = 0; i < np; i++) { + for (int j = 0; j < nps + 2; j++) { + fscanf(f, "%lg", &pr[j]); + fprintf(g, "%.10le ", pr[j]); + } + fprintf(g, "\n"); + } + fclose(f); + + scanbumper = bumperlist; + for (il = 1; il <= nmod; il++) { + for (int i = 0; i < nps; i++) { + pr[i] = scanbumper->p0[i]; + } + double u0 = pr[2]; + double alpha = pr[3]; + double s0, s = exp(pr[0]), s2; + double q = exp(pr[1]); + double q2 = 0.000001; + /*double dt = (scanbumper->tanomaly - pr[6]) / exp(pr[5]);*/ + /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ + double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); + double xc0, xc; + double rho = exp(pr[4]); + double y1A, y2A; + double beta, beta0; + + //y1A e y2A posizione a tanomaly + //mentre y1anomaly e y2anomaly posizione ricavate dal LevMarfit + y1A = pr[2] * sin(pr[3]) - dt * cos(pr[3]) + q * s / (1 + q); //rispetto la primaria + y2A = -pr[2] * cos(pr[3]) - dt * sin(pr[3]); + beta = beta0 = atan2(y2A, y1A); + + xc0 = sqrt(y1A * y1A + y2A * y2A); + xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) + xc); + while (xc < 4 * sqrt(q2) / (s0 * s0)) q2 *= 0.1; + xc = xc0 + 4 * sqrt(q2) / (s0 * s0); + s2 = 0.5 * (sqrt(4 + xc * xc) + xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + xc = xc0 - 4 * sqrt(q2) / (s0 * s0); + s2 = 0.5 * (sqrt(4 + xc * xc) + xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + + xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) - xc); + q2 = 0.00001; + while (xc < 3 * sqrt(3 * q2) * s0 * s0 * s0) q2 *= 0.1; + xc = xc0 - 3 * sqrt(3 * q2) * s0 * s0 * s0; + s2 = 0.5 * (sqrt(4 + xc * xc) - xc); + beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + xc = xc0 + 3 * sqrt(3 * q2) * s0 * s0 * s0; + s2 = 0.5 * (sqrt(4 + xc * xc) - xc); + beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + + scanbumper = scanbumper->next; + } + fclose(g); + free(peaks); + remove("InitCondTX.txt"); + rename("InitCondTX-temp.txt", "InitCondTX.txt"); + } + } + if (modelcode[1] == 'O') { printf("\n- Preparing initial conditions for Keplerian orbital motion"); if (f = fopen("InitCondLK.txt", "r")) { @@ -805,6 +1017,7 @@ int main(int argc, char* argv[]) { rename("InitCondLK-temp.txt", "InitCondLK.txt"); } } + break; case 'P': if (modelcode[1] == 'S') { @@ -892,8 +1105,9 @@ int main(int argc, char* argv[]) { /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { if (ipeak == idpeak) continue;*/ double q = 0.001; + /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); + double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3]));*/ double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); - /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ double xc0 = sqrt(u0 * u0 + dt * dt), xc; double alpha0 = atan2(u0, -dt), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); @@ -993,10 +1207,11 @@ int main(int argc, char* argv[]) { /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { if (ipeak == idpeak) continue;*/ double q = 0.001; - double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); - //double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); + /*double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]);*/ + double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); double xc0 = sqrt(u0 * u0 + dt * dt), xc; - double alpha0 = atan2(u0, -dt), alpha; + //double alpha0 = atan2(u0, -dt), alpha; + double alpha0 = atan2(scanbumper->y2anomaly-u0, scanbumper->y1anomaly), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); @@ -1108,25 +1323,16 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } - /*double mindpeak = 1.e100; - int idpeak = 0; - for (int ipeak = 0; ipeak < npeaks; ipeak++) { - double ddpeak = fabs(peaks[ipeak] - pr[2]); - if (ddpeak < mindpeak) { - idpeak = ipeak; - mindpeak = ddpeak; - } - }*/ - double u0 = pr[0]; double s0, s; /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { if (ipeak == idpeak) continue;*/ double q = 0.001; - double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); - //double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); + //double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); + double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); double xc0 = sqrt(u0 * u0 + dt * dt), xc; - double alpha0 = atan2(u0, -dt), alpha; + double alpha0 = atan2(scanbumper->y2anomaly - u0, scanbumper->y1anomaly), alpha; + /*double alpha0 = atan2(u0, -dt); */ double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); //funzione inline @@ -1190,6 +1396,7 @@ int main(int argc, char* argv[]) { } } break; + case 'B': if (modelcode[1] == 'S') { // Preparation of initial conditions for xallarap From 589b7e99693dae79a3b3e225d8f642b1f54bde7b Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Tue, 20 Jan 2026 13:17:24 +0100 Subject: [PATCH 25/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 215 +++++++++++++++++----------------- 1 file changed, 106 insertions(+), 109 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index 9616b48..6312563 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -726,107 +726,104 @@ int main(int argc, char* argv[]) { } } } - } - printf("\n- Adding initial conditions for planets"); - if (f = fopen("InitCondTS.txt", "r")) { - int npeaks = 0; - double* peaks; - g = fopen("InitCondTS-temp.txt", "w"); - fscanf(f, "%d %d", &npeaks, &np); - if (npeaks == 0) { + + printf("\n- Adding initial conditions for planets"); + if (f = fopen("InitCondTS.txt", "r")) { + int npeaks = 0; + double* peaks; + g = fopen("InitCondTS-temp.txt", "w"); + fscanf(f, "%d %d", &npeaks, &np); + if (npeaks == 0) { fclose(g); fclose(f); remove("InitCondTS-temp.txt"); - } - else { - fprintf(g, "%d %d\n", npeaks, np + ((npeaks > 0) ? 6 * nmod : 0)); - peaks = (double*)malloc(sizeof(double) * npeaks); - - printf("\nNumber of initial conditions: %d", np + ((npeaks > 0) ? 6 * nmod : 0)); - for (int i = 0; i < npeaks; i++) { - fscanf(f, "%lg", &peaks[i]); - fprintf(g, "%le", peaks[i]); - fscanf(f, "%[^\n]s", &buffer); - fprintf(g, "%s\n", buffer); } - for (int i = 0; i < np; i++) { - for (int j = 0; j < 10; j++) { - fscanf(f, "%lg", &pr[j]); - fprintf(g, "%.10le ", pr[j]); - } - fprintf(g, "\n"); - } - fclose(f); - + else { + fprintf(g, "%d %d\n", npeaks, np + ((npeaks > 0) ? 6 * nmod : 0)); + peaks = (double*)malloc(sizeof(double) * npeaks); - scanbumper = bumperlist; - for (il = 1; il <= nmod; il++) { - for (int i = 0; i < nps; i++) { - pr[i] = scanbumper->p0[i]; + printf("\nNumber of initial conditions: %d", np + ((npeaks > 0) ? 6 * nmod : 0)); + for (int i = 0; i < npeaks; i++) { + fscanf(f, "%lg", &peaks[i]); + fprintf(g, "%le", peaks[i]); + fscanf(f, "%[^\n]s", &buffer); + fprintf(g, "%s\n", buffer); } + for (int i = 0; i < np; i++) { + for (int j = 0; j < 10; j++) { + fscanf(f, "%lg", &pr[j]); + fprintf(g, "%.10le ", pr[j]); + } + fprintf(g, "\n"); + } + fclose(f); - double u0 = pr[2]; - double alpha = pr[3]; - double s0, s = exp(pr[0]), s2; - double q = exp(pr[1]); - double q2 = 0.000001; - /*double dt = (scanbumper->tanomaly - pr[6]) / exp(pr[5]);*/ - /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ - double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); - double xc0, xc; - double rho = exp(pr[4]); - double y1A, y2A; - double beta, beta0; - - //y1A e y2A posizione a tanomaly - //mentre y1anomaly e y2anomaly posizione ricavate dal LevMarfit - y1A = pr[2] * sin(pr[3]) - dt * cos(pr[3]) + q * s / (1 + q); //rispetto la primaria - y2A = -pr[2] * cos(pr[3]) - dt * sin(pr[3]); - beta = beta0 = atan2(y2A, y1A); - - xc0 = sqrt(y1A * y1A + y2A * y2A); - xc = xc0; - s0 = 0.5 * (sqrt(4 + xc * xc) + xc); - while (xc < 4 * sqrt(q2) / (s0 * s0)) q2 *= 0.1; - xc = xc0 + 4 * sqrt(q2) / (s0 * s0); - s2 = 0.5 * (sqrt(4 + xc * xc) + xc); - - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); - xc = xc0 - 4 * sqrt(q2) / (s0 * s0); - s2 = 0.5 * (sqrt(4 + xc * xc) + xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); - - xc = xc0; - s0 = 0.5 * (sqrt(4 + xc * xc) - xc); - q2 = 0.00001; - while (xc < 3 * sqrt(3 * q2) * s0 * s0 * s0) q2 *= 0.1; - xc = xc0 - 3 * sqrt(3 * q2) * s0 * s0 * s0; - s2 = 0.5 * (sqrt(4 + xc * xc) - xc); - beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); - - - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); - beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); - - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); - xc = xc0 + 3 * sqrt(3 * q2) * s0 * s0 * s0; - s2 = 0.5 * (sqrt(4 + xc * xc) - xc); - beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); - beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + scanbumper = bumperlist; + for (il = 1; il <= nmod; il++) { + for (int i = 0; i < nps; i++) { + pr[i] = scanbumper->p0[i]; + } - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + double u0 = pr[2]; + double alpha = pr[3]; + double s0, s = exp(pr[0]), s2; + double q = exp(pr[1]); + double q2 = 0.000001; + /*double dt = (scanbumper->tanomaly - pr[6]) / exp(pr[5]);*/ + /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ + double xc0, xc; + double rho = exp(pr[4]); + double y1A = scanbumper->y1anomaly + q * s / (1 + q), y2A = scanbumper->y2anomaly; + double beta, beta0; + + //y1A = pr[2] * sin(pr[3]) - dt * cos(pr[3]) + q * s / (1 + q); //rispetto la primaria + //y2A = -pr[2] * cos(pr[3]) - dt * sin(pr[3]); + beta = beta0 = atan2(y2A, y1A); + + xc0 = sqrt(y1A * y1A + y2A * y2A); + xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) + xc); + while (xc < 4 * sqrt(q2) / (s0 * s0)) q2 *= 0.1; + xc = xc0 + 4 * sqrt(q2) / (s0 * s0); + s2 = 0.5 * (sqrt(4 + xc * xc) + xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + xc = xc0 - 4 * sqrt(q2) / (s0 * s0); + s2 = 0.5 * (sqrt(4 + xc * xc) + xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + + xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) - xc); + q2 = 0.00001; + while (xc < 3 * sqrt(3 * q2) * s0 * s0 * s0) q2 *= 0.1; + xc = xc0 - 3 * sqrt(3 * q2) * s0 * s0 * s0; + s2 = 0.5 * (sqrt(4 + xc * xc) - xc); + beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + xc = xc0 + 3 * sqrt(3 * q2) * s0 * s0 * s0; + s2 = 0.5 * (sqrt(4 + xc * xc) - xc); + beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); + beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); - scanbumper = scanbumper->next; + scanbumper = scanbumper->next; + } + fclose(g); + free(peaks); + remove("InitCondTS.txt"); + rename("InitCondTS-temp.txt", "InitCondTS.txt"); } - fclose(g); - free(peaks); - remove("InitCondTS.txt"); - rename("InitCondTS-temp.txt", "InitCondTS.txt"); } } - if (modelcode[1] == 'X') { printf("\n- Preparing initial conditions for orbital motion"); if (f = fopen("InitCondLO.txt", "r")) { @@ -890,8 +887,6 @@ int main(int argc, char* argv[]) { rename("InitCondLK-temp.txt", "InitCondLK.txt"); } } - - ////bisogna aggiungere TX printf("\n- Adding initial conditions for planets"); if (f = fopen("InitCondTX.txt", "r")) { @@ -930,16 +925,13 @@ int main(int argc, char* argv[]) { double q2 = 0.000001; /*double dt = (scanbumper->tanomaly - pr[6]) / exp(pr[5]);*/ /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ - double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); double xc0, xc; double rho = exp(pr[4]); - double y1A, y2A; + double y1A = scanbumper->y1anomaly + q * s / (1 + q), y2A = scanbumper->y2anomaly; double beta, beta0; - //y1A e y2A posizione a tanomaly - //mentre y1anomaly e y2anomaly posizione ricavate dal LevMarfit - y1A = pr[2] * sin(pr[3]) - dt * cos(pr[3]) + q * s / (1 + q); //rispetto la primaria - y2A = -pr[2] * cos(pr[3]) - dt * sin(pr[3]); + //y1A = pr[2] * sin(pr[3]) - dt * cos(pr[3]) + q * s / (1 + q); //rispetto la primaria + //y2A = -pr[2] * cos(pr[3]) - dt * sin(pr[3]); beta = beta0 = atan2(y2A, y1A); xc0 = sqrt(y1A * y1A + y2A * y2A); @@ -1102,18 +1094,21 @@ int main(int argc, char* argv[]) { double u0 = exp(pr[0]); double s0, s; + double y1A = scanbumper->y1anomaly, y2A = scanbumper->y2anomaly; /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { if (ipeak == idpeak) continue;*/ double q = 0.001; /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]); - double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3]));*/ - double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); - double xc0 = sqrt(u0 * u0 + dt * dt), xc; - double alpha0 = atan2(u0, -dt), alpha; + double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]);*/ + //double xc0 = sqrt(u0 * u0 + dt * dt), xc; + double xc0 = sqrt(y1A * y1A + y2A * y2A), xc; + //double alpha0 = atan2(u0, -dt), alpha; + double alpha0 = atan2(-y2A, y1A), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); xc = xc0; + s0 = 0.5 * (sqrt(4 + xc * xc) + xc); while (xc < 4 * sqrt(q) / (s0 * s0)) q *= 0.1; xc = xc0 + 4 * sqrt(q) / (s0 * s0); @@ -1201,17 +1196,18 @@ int main(int argc, char* argv[]) { mindpeak = ddpeak; } }*/ - + double u0 = pr[0]; double s0, s; + double y1A = scanbumper->y1anomaly, y2A = scanbumper->y2anomaly; /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { if (ipeak == idpeak) continue;*/ double q = 0.001; /*double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]);*/ - double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); - double xc0 = sqrt(u0 * u0 + dt * dt), xc; + double xc0 = sqrt(y1A * y1A + y2A * y2A), xc; + double alpha0 = atan2(-y2A, y1A), alpha; + //double xc0 = sqrt(u0 * u0 + dt * dt), xc; //double alpha0 = atan2(u0, -dt), alpha; - double alpha0 = atan2(scanbumper->y2anomaly-u0, scanbumper->y1anomaly), alpha; double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); @@ -1322,17 +1318,18 @@ int main(int argc, char* argv[]) { for (int i = 0; i < nps; i++) { pr[i] = scanbumper->p0[i]; } - + double u0 = pr[0]; double s0, s; + double y1A = scanbumper->y1anomaly, y2A = scanbumper->y2anomaly; /*for (int ipeak = 0; ipeak < npeaks; ipeak++) { if (ipeak == idpeak) continue;*/ double q = 0.001; //double dt = (scanbumper->tanomaly - pr[2]) / exp(pr[1]); - double dt = -(scanbumper->y2anomaly + pr[2] * cos(pr[3]) / sin(pr[3])); - double xc0 = sqrt(u0 * u0 + dt * dt), xc; - double alpha0 = atan2(scanbumper->y2anomaly - u0, scanbumper->y1anomaly), alpha; - /*double alpha0 = atan2(u0, -dt); */ + double xc0 = sqrt(y1A * y1A + y2A * y2A), xc; + double alpha0 = atan2(-y2A, y1A), alpha; + /*double xc0 = sqrt(u0 * u0 + dt * dt), xc; + double alpha0 = atan2(u0, -dt), alpha; */ double rho = 0.001; // exp(pr[3] - sqrt(scanbumper->cov[3 * nps + 3])); //funzione inline From bb193e938ac922e42724b189c5baff2429f4c8ba Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Tue, 20 Jan 2026 13:18:15 +0100 Subject: [PATCH 26/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 65 +++++++++++++++++++++++---------------- 1 file changed, 38 insertions(+), 27 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 1595b83..3316ed9 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -1,6 +1,6 @@ // LevMarFit.cpp // Implementation of all methods in LevMarFit.h for Levenberg-Marquardt fitting - +// ACversion #define _CRT_SECURE_NO_WARNINGS #define _USE_MATH_DEFINES #include "LevMarFit.h" @@ -14,6 +14,7 @@ #include #include + //using namespace std; using std::regex, std::string, std::regex_match; using namespace std::filesystem; @@ -55,6 +56,7 @@ std::vector> logposs = { {0, 1, 3}, const double epsilon = 1.e-100; + LevMar::LevMar(int argc, char* argv[]) { setbuf(stdout, nullptr); printf("******************************************\n"); @@ -72,6 +74,7 @@ LevMar::LevMar(int argc, char* argv[]) { VBM->RelTol = 0.001; VBM->parallaxsystem = 1; VBM->SetMethod(VBMicrolensing::Method::Multipoly); + ReadFiles(argc, argv); @@ -188,8 +191,6 @@ void LevMar::ReadFiles(int argc, char* argv[]) { printf("\n\n- Event: %s\n", eventname); current_path(eventname); - printf("- Working directory: %s\n", current_path().string().c_str()); - //if(argc>2){ // strcpy(eventname,argv[1]); // strcpy(outdir,argv[2]); @@ -215,6 +216,7 @@ void LevMar::ReadFiles(int argc, char* argv[]) { break; } } + current_path(eventname); printf("\n- Event directory: %s\n", current_path().string().c_str()); @@ -238,14 +240,14 @@ void LevMar::ReadFiles(int argc, char* argv[]) { error = InitCond(presigmapr, preleftlim, prerightlim); pr[1] = log(pr[1]); pr[3] = log(pr[3]); - //current_path(exedir); - //current_path(".."); + /*current_path(exedir);*/ + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); - - } + getchar(); + } else { modnumber = 0; nps = 4; @@ -261,13 +263,10 @@ void LevMar::ReadFiles(int argc, char* argv[]) { } //current_path(exedir); //current_path(".."); - current_path("data"); VBM->LoadESPLTable("ESPL.tbl"); current_path(eventname); - break; - case 'B': if (modelcode[1] == 'O') { modnumber = 3; @@ -282,7 +281,6 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[0] = log(pr[0]); pr[1] = log(pr[1]); pr[6] = log(pr[6]); - //current_path(exedir); //current_path(".."); current_path("data"); @@ -407,7 +405,6 @@ void LevMar::ReadFiles(int argc, char* argv[]) { current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); - } else { modnumber = 8; @@ -428,7 +425,6 @@ void LevMar::ReadFiles(int argc, char* argv[]) { break; } printf("\n- Model: %s", modelcode); - } catch (...) { error = 10; @@ -693,7 +689,6 @@ void LevMar::ReadOptions(double* preleftlim, double* prerightlim, double* presig VBM->lens_mass_luminosity_exponent = lens_mass_luminosity_exponent; current_path(eventname); } - int LevMar::InitCond(double* presigmapr, double* preleftlim, double* prerightlim) { char buffer[3200], initcondfile[256]; int npeaks, ninit, incond; @@ -703,6 +698,7 @@ int LevMar::InitCond(double* presigmapr, double* preleftlim, double* prerightlim rightlim = (double*)malloc(sizeof(double) * nps); pr = (double*)malloc(sizeof(double) * nps); if (parametersfile[0] == 0) { + sprintf(initcondfile, "InitCond%c%c.txt", modelcode[0], modelcode[1]); current_path("InitCond"); f = fopen(initcondfile, "r"); @@ -781,7 +777,7 @@ void LevMar::ReadAncillary() { } current_path("Data"); - + // If Normalization.txt is present, use numbers therein to normalize datasets. normfacs = (double*)malloc(sizeof(double) * nfil); @@ -1216,6 +1212,7 @@ int LevMar::Run() { // If this fit has found the best minimum so far, update minchi.dat // current_path(".."); + strcpy(filename, "minchi.dat"); strcpy(filename, "minchi.dat"); if (exists(filename)) { f = fopen(filename, "r"); @@ -1294,9 +1291,6 @@ void LevMar::EvaluateModel(double* pr, int fl, int ips) { VBM->BinaryLightCurveKepler(pr, tfl, fbfl, y1a, y2a, seps, sizes[fl]); } break; - - - case 8: VBM->TripleLightCurve(pr, tfl, fbfl, y1a, y2a, sizes[fl]); break; @@ -1309,7 +1303,6 @@ void LevMar::EvaluateModel(double* pr, int fl, int ips) { } break; } - } double LevMar::ChiSquared(double* pr) { @@ -1317,7 +1310,7 @@ double LevMar::ChiSquared(double* pr) { double p1max = 0, maxsump = 0; double maxsumn = 0; double t1 = 0, t2 = 0, tmax = 0; - + double y1max=0, y2max=0, y1t1=0, y2t1=0, y1t2=0, y2t2=0; maxmaxsum = 0; for (int fl = 0; fl < nfil; fl++) { @@ -1418,6 +1411,9 @@ double LevMar::ChiSquared(double* pr) { if (p1 > p1max) { //per trovare il picco p1max = p1; // aggiorna il massimo residuo positivo tmax = t[i]; // tempo in cui si ha il massimo residuo positivo + y1max = y1a[i]; + y2max = y2a[i]; + } } @@ -1426,10 +1422,13 @@ double LevMar::ChiSquared(double* pr) { if (in_pos_sequence) { in_pos_sequence = false; t1 = t[i-1]; // tempo ultimo pos prima di una sequenza neg + y1t1 = y1a[i - 1]; + y2t1 = y2a[i - 1]; } ////aggiorna t2 t2 = t[i]; - + y1t2 = y1a[i]; + y2t2 = y2a[i]; maxsump = 0; // azzera la somma dei residui positivi maxsumn += p1; // somma dei residui negativi, è un numero negativo p1max = 0; @@ -1439,16 +1438,21 @@ double LevMar::ChiSquared(double* pr) { if (maxsump > maxmaxsum) { maxmaxsum = maxsump; tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi - + y1maxmax = y1max; + y2maxmax = y2max; } else if(fabs(maxsumn) > maxmaxsum); { maxmaxsum = - maxsumn; //calcolo del tempo || pr[6] è t0 in binary lens if (abs(double(pr[6] - t1)) > abs(double(pr[6] - t2))) { tmaxmax = t1; + y1maxmax = y1t1; + y2maxmax = y2t1; } else { tmaxmax = t2; + y1maxmax = y1t2; + y2maxmax = y2t2; } } } @@ -1463,15 +1467,20 @@ double LevMar::ChiSquared(double* pr) { if (maxsump > maxmaxsum) { maxmaxsum = maxsump; tmaxmax = tmax; // tempo in cui si ha la somma massima dei residui positivi consecutivi - + y1maxmax = y1max; + y2maxmax = y2max; } else if (fabs(maxsumn) > maxmaxsum); { maxmaxsum = - maxsumn; if (abs(pr[6] - t1) > abs(pr[6] - t2)) { tmaxmax = t1; + y1maxmax = y1t1; + y2maxmax = y2t1; } else { tmaxmax = t2; + y1maxmax = y1t2; + y2maxmax = y2t2; } } @@ -1670,8 +1679,6 @@ void LevMar::Covariance() { } } - - void LevMar::PrintOut(double* pr) { int npp = nps, ilog = 0, logsize = logposs[modnumber].size(); double fl; @@ -1847,8 +1854,12 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { //stampa tmaxmax fprintf(f, "%.16le ", tmaxmax); - //stampa maxmaxsum - fprintf(f, "%.16le ", maxmaxsum); + //stampa y1maxmax e y2maxmax + fprintf(f, "%.16le\n", y1maxmax); + fprintf(f, "%.16le\n", y2maxmax); + + //stampa delta chi quadro + fprintf(f, "%.16le", maxmaxsum); // Write chi square fprintf(f, "%.16le\n", c0); From 68716c59aeab0d23dd517880bb59fe26f7b1ef0e Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Tue, 20 Jan 2026 13:19:01 +0100 Subject: [PATCH 27/40] Add files via upload --- RTModel/include/bumper.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RTModel/include/bumper.h b/RTModel/include/bumper.h index 295f5f5..8c4b167 100644 --- a/RTModel/include/bumper.h +++ b/RTModel/include/bumper.h @@ -14,6 +14,8 @@ class bumper{ double *curv, *cov; double Amp; double tanomaly; + double y1anomaly, y2anomaly; + double y1maxmax, y2maxmax; double resanomaly; double maxsum, maxmaxsum, tmaxmax, tmax; int nps; From 220195538268cb5552c28b5837468c2588bf4ba5 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 28 Jan 2026 10:05:42 +0100 Subject: [PATCH 28/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 58 ++++++++++++++++++++------------------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 3316ed9..82812e5 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -240,13 +240,12 @@ void LevMar::ReadFiles(int argc, char* argv[]) { error = InitCond(presigmapr, preleftlim, prerightlim); pr[1] = log(pr[1]); pr[3] = log(pr[3]); - /*current_path(exedir);*/ + current_path(exedir); current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); - getchar(); } else { modnumber = 0; @@ -261,11 +260,12 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[3] = log(pr[3]); } - //current_path(exedir); - //current_path(".."); - current_path("data"); - VBM->LoadESPLTable("ESPL.tbl"); - current_path(eventname); + current_path(exedir); + current_path(".."); + current_path("data"); + VBM->LoadESPLTable("ESPL.tbl"); + current_path(eventname); + break; case 'B': if (modelcode[1] == 'O') { @@ -281,8 +281,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[0] = log(pr[0]); pr[1] = log(pr[1]); pr[6] = log(pr[6]); - //current_path(exedir); - //current_path(".."); + current_path(exedir); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -300,12 +300,12 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[0] = log(pr[0]); pr[1] = log(pr[1]); pr[6] = log(pr[6]); - } - //current_path(exedir); - //current_path(".."); - current_path("data"); - VBM->LoadESPLTable("ESPL.tbl"); - current_path(eventname); + } + current_path(exedir); + current_path(".."); + current_path("data"); + VBM->LoadESPLTable("ESPL.tbl"); + current_path(eventname); break; case 'L': if (modelcode[1] == 'X') { @@ -321,8 +321,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[4] = log(pr[4]); pr[5] = log(pr[5]); - //current_path(exedir); - //current_path(".."); + /*current_path(exedir); + current_path("..");*/ current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -341,8 +341,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[4] = log(pr[4]); pr[5] = log(pr[5]); - //current_path(exedir); - //current_path(".."); + current_path(exedir); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -361,8 +361,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[4] = log(pr[4]); pr[5] = log(pr[5]); - //current_path(exedir); - //current_path(".."); + current_path(exedir); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -400,8 +400,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[5] = log(pr[5]); pr[7] = log(pr[7]); pr[8] = log(pr[8]); - //current_path(exedir); - //current_path(".."); + current_path(exedir); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -584,7 +584,7 @@ void LevMar::ReadOptions(double* preleftlim, double* prerightlim, double* presig } if (command[0] == 'g' && command[1] == '_') { // Blending - sscanf(&command[2], "%d", &consindex[conscurrent]); + scanf(&command[2], "%d", &consindex[conscurrent]); consindex[conscurrent] += 10000; } // special combinations @@ -1855,11 +1855,11 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { fprintf(f, "%.16le ", tmaxmax); //stampa y1maxmax e y2maxmax - fprintf(f, "%.16le\n", y1maxmax); - fprintf(f, "%.16le\n", y2maxmax); + fprintf(f, "%.16le ", y1maxmax); + fprintf(f, "%.16le ", y2maxmax); //stampa delta chi quadro - fprintf(f, "%.16le", maxmaxsum); + fprintf(f, "%.16le ", maxmaxsum); // Write chi square fprintf(f, "%.16le\n", c0); @@ -1889,7 +1889,7 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { fprintf(f, " %le", errs[i]); } fprintf(f, "\n"); - + // Print covariance matrix to file for (int i = 0; i < nps; i++) { if (i > 0) fprintf(f, "\n"); @@ -1899,6 +1899,8 @@ void LevMar::PrintFile(char* filename, int il, double c0, bool printerrors) { } } fprintf(f, "\n"); + } + fclose(f); } \ No newline at end of file From 41b435e5d85e4d58d7286c1b0e3abacb937da79e Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Mon, 23 Feb 2026 12:13:09 +0100 Subject: [PATCH 29/40] Add files via upload --- RTModel/lib/InitCond.cpp | 180 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 177 insertions(+), 3 deletions(-) diff --git a/RTModel/lib/InitCond.cpp b/RTModel/lib/InitCond.cpp index 5e67ade..560d17f 100644 --- a/RTModel/lib/InitCond.cpp +++ b/RTModel/lib/InitCond.cpp @@ -30,7 +30,7 @@ bool onlyorbital = false; // Only orbital motion models will be calculated. bool onlyupdate = false; // No model search, but only update of previously found best models. int usesatellite = 0; // Satellite to be used for initial conditions. Ground telescopes by default. char templatelibrary[256] = ""; // User-specified template library -char modelcategories[256] = "PSPXBSBOLSLXLO"; +char modelcategories[256] = "PSPXBSBOLSLXLOTXTS"; char astroini[256] = "0.0 0.0 0.125 1.0"; //double tau = 0.5; // Provisional!!! Exclude peaks shorter than tau @@ -239,6 +239,11 @@ int main(int argc, char* argv[]) pfol = prm + 2; strcpy(prm, pfol); } + prm = strstr(modelcategories, "TS"); + if (prm != 0) { + pfol = prm + 2; + strcpy(prm, pfol); + } if (onlyorbital) { prm = strstr(modelcategories, "LX"); if (prm != 0) { @@ -246,6 +251,13 @@ int main(int argc, char* argv[]) strcpy(prm, pfol); } } + if (onlyorbital) { + prm = strstr(modelcategories, "TX"); + if (prm != 0) { + pfol = prm + 2; + strcpy(prm, pfol); + } + } } // If only update of previous models is required, there is no need to calculate peaks of current data @@ -894,7 +906,8 @@ int main(int argc, char* argv[]) else maxoldmodels = 0; current_path(eventname); - if (!exists("InitCond")) create_directory("InitCond"); + if (!exists("InitCond")) + create_directory("InitCond"); current_path("InitCond"); searchstring = regex(".*Init.*\\.txt"); @@ -1416,7 +1429,7 @@ int main(int argc, char* argv[]) fclose(f); free(tt); current_path(eventname); - } + } dn = 0; if (strstr(modelcategories, "LO") != 0) { @@ -1652,6 +1665,167 @@ int main(int argc, char* argv[]) free(tt); } + dn = 0; + if (strstr(modelcategories, "TS") != 0) { + filebest = regex("TS.*\\.txt"); + strcpy(fileinit, "InitCondTS.txt"); + nps = npsold = 10; + if (astrometric) nps += 4; + if (astrometricold && astrometric) npsold += 4; + + tt = (double*)malloc(sizeof(double) * nps * maxoldmodels); + if (exists(path(runstring) / path("Models"))) { + current_path(path(runstring) / path("Models")); + for (auto const& itr : directory_iterator(".")) { + if (dn >= maxoldmodels) break; + string curfile = (itr).path().filename().string(); + if (regex_match(curfile, filebest)) { + f = fopen(curfile.c_str(), "r"); + for (int j = 0; j < npsold; j++) { + if (fscanf(f, "%le", &tt[dn * nps + j]) < 1) { + sscanf(astroini, "%lf %lf %lf %lf", &(tt[dn * nps + j]), &(tt[dn * nps + j + 1]), &(tt[dn * nps + j + 2]), &(tt[dn * nps + j + 3])); + break; + } + } + fclose(f); + dn++; + } + } + } + printf("\n- Writing initial conditions for fitting to %s\n\n", fileinit); + current_path(eventname); + current_path("InitCond"); + + f = fopen(fileinit, "w"); + fprintf(f, "0 %d\n", dn); + + // First we write the initial conditions from previous best models + for (int i = 0; i < dn; i++) { + for (int j = 0; j < nps; j++) { + fprintf(f, "%le ", tt[i * nps + j]); + } + if (npsold < nps) fprintf(f, "%s", astroini); + fprintf(f, "\n"); + } + fclose(f); + free(tt); + current_path(eventname); + } + + dn = 0; + + if (strstr(modelcategories, "TX") != 0) { + filebest = regex("TX.*\\.txt"); + strcpy(fileinit, "InitCondTX.txt"); + nps = npsold = 11; + if (astrometric) nps += 4; + if (astrometricold && astrometric) npsold += 4; + + tt = (double*)malloc(sizeof(double) * nps * maxoldmodels); + if (exists(path(runstring) / path("Models"))) { + current_path(path(runstring) / path("Models")); + for (auto const& itr : directory_iterator(".")) { + if (dn >= maxoldmodels) break; + string curfile = (itr).path().filename().string(); + if (regex_match(curfile, filebest)) { + f = fopen(curfile.c_str(), "r"); + for (int j = 0; j < npsold; j++) { + if (fscanf(f, "%le", &tt[dn * nps + j]) < 1) { + sscanf(astroini, "%lf %lf %lf %lf", &(tt[dn * nps + j]), &(tt[dn * nps + j + 1]), &(tt[dn * nps + j + 2]), &(tt[dn * nps + j + 3])); + break; + } + } + fclose(f); + dn++; + } + } + } + + printf("\n- Writing initial conditions for fitting to %s\n\n", fileinit); + current_path(eventname); + current_path("InitCond"); + + f = fopen(fileinit, "w"); + if (strstr(modelcategories, "TS") != 0) { + fprintf(f, "0 %d\n", dn); + } + //else { + // fprintf(f, "%d %d\n", newpeaks->length, np * (newpeaks->length * (newpeaks->length - 1)) * 2 + dn); + // // Then we write the characteristics of the peaks used + // for (p = newpeaks->first; p; p = p->next) { + // fprintf(f, "%le %le %le %le %le\n", p->t, p->tl, p->tr, p->y, p->sig); + // } + //} + + // First we write the initial conditions from previous best models + for (int i = 0; i < dn; i++) { + for (int j = 0; j < nps; j++) { + fprintf(f, "%le ", tt[i * nps + j]); + } + if (npsold < nps) fprintf(f, "%s", astroini); + fprintf(f, "\n"); + } + // Here we write the initial conditions by matching the newpeaks to the peaks recorded in the template library + //if (strstr(modelcategories, "TS") == 0) { + // if (newpeaks->length > 1) { + // for (int i = 0; i < np; i++) { + // for (pl = newpeaks->first; pl->next; pl = pl->next) { + // for (pr = pl->next; pr; pr = pr->next) { + // t1 = (pl->t < pr->t) ? pl->t : pr->t; + // t2 = (pl->t < pr->t) ? pr->t : pl->t; + // tE = (t2 - t1) / (yy[i * 7 + 6] - yy[i * 7 + 5]); // yy[i*7+6] and yy[i*7+5] are the peak times of the ith template + // t0 = t2 - tE * yy[i * 7 + 6]; + // for (int j = 0; j < 5; j++) { + // fprintf(f, "%le ", yy[i * 7 + j]); // We use the s,q,u0,alpha,rho parameters from the template + // } + // fprintf(f, "%le %le", tE, t0); // and use tE and t0 from the time matching + // fprintf(f, " 0.0 0.0"); // parallax for nostatic + // if (astrometric) fprintf(f, " %s", astroini); + // fprintf(f, "\n"); + + // // Reflected initial condition for nostatic + // yy[i * 7 + 2] = -yy[i * 7 + 2]; + // yy[i * 7 + 3] = -yy[i * 7 + 3]; + // for (int j = 0; j < 5; j++) { + // fprintf(f, "%le ", yy[i * 7 + j]); + // } + // fprintf(f, "%le %le", tE, t0); + // fprintf(f, " 0.0 0.0"); // parallax for nostatic + // if (astrometric) fprintf(f, " %s", astroini); + // fprintf(f, "\n"); + + // tE = (t1 - t2) / (yy[i * 7 + 6] - yy[i * 7 + 5]); // We also include the time-reverse matching + // t0 = t1 - tE * yy[i * 7 + 6]; + // yy[i * 7 + 2] = -yy[i * 7 + 2]; //u0 and alpha are reversed + // yy[i * 7 + 3] = yy[i * 7 + 3] + M_PI; + // tE = -tE; + // for (int j = 0; j < 5; j++) { + // fprintf(f, "%le ", yy[i * 7 + j]); + // } + // fprintf(f, "%le %le", tE, t0); // and use tE and t0 from the time matching + // fprintf(f, " 0.0 0.0"); // parallax for nostatic + // if (astrometric) fprintf(f, " %s", astroini); + // fprintf(f, "\n"); + + // // Reflected initial condition for nostatic + // yy[i * 7 + 2] = -yy[i * 7 + 2]; + // yy[i * 7 + 3] = -yy[i * 7 + 3]; + // for (int j = 0; j < 5; j++) { + // fprintf(f, "%le ", yy[i * 7 + j]); + // } + // fprintf(f, "%le %le", tE, t0); + // fprintf(f, " 0.0 0.0"); // parallax for nostatic + // if (astrometric) fprintf(f, " %s", astroini); + // fprintf(f, "\n"); + // } + // } + // } + // } + //} + fclose(f); + free(tt); + current_path(eventname); + } printf("\n---- Done"); //Sleep(5000l); From c088232f8602e066ac2fdc489651af488ccec6ac Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Mon, 23 Feb 2026 14:29:19 +0100 Subject: [PATCH 30/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index 6312563..713e487 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -733,16 +733,16 @@ int main(int argc, char* argv[]) { double* peaks; g = fopen("InitCondTS-temp.txt", "w"); fscanf(f, "%d %d", &npeaks, &np); - if (npeaks == 0) { + /*if (npeaks == 0) { fclose(g); fclose(f); remove("InitCondTS-temp.txt"); - } - else { - fprintf(g, "%d %d\n", npeaks, np + ((npeaks > 0) ? 6 * nmod : 0)); - peaks = (double*)malloc(sizeof(double) * npeaks); + }*/ + /*else {*/ + fprintf(g, "%d %d\n", npeaks, np + 6 * nmod); + peaks = (double*)malloc(sizeof(double) * (npeaks + 1)); - printf("\nNumber of initial conditions: %d", np + ((npeaks > 0) ? 6 * nmod : 0)); + printf("\nNumber of initial conditions: %d", np + 6 * nmod ); for (int i = 0; i < npeaks; i++) { fscanf(f, "%lg", &peaks[i]); fprintf(g, "%le", peaks[i]); @@ -769,7 +769,7 @@ int main(int argc, char* argv[]) { double alpha = pr[3]; double s0, s = exp(pr[0]), s2; double q = exp(pr[1]); - double q2 = 0.000001; + double q2 = 0.001; /*double dt = (scanbumper->tanomaly - pr[6]) / exp(pr[5]);*/ /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ double xc0, xc; @@ -821,7 +821,7 @@ int main(int argc, char* argv[]) { free(peaks); remove("InitCondTS.txt"); rename("InitCondTS-temp.txt", "InitCondTS.txt"); - } + //} } } if (modelcode[1] == 'X') { @@ -922,7 +922,7 @@ int main(int argc, char* argv[]) { double alpha = pr[3]; double s0, s = exp(pr[0]), s2; double q = exp(pr[1]); - double q2 = 0.000001; + double q2 = 0.001; /*double dt = (scanbumper->tanomaly - pr[6]) / exp(pr[5]);*/ /*double dt = (peaks[ipeak] - pr[2]) / exp(pr[1]);*/ double xc0, xc; @@ -1050,16 +1050,16 @@ int main(int argc, char* argv[]) { double* peaks; g = fopen("InitCondLS-temp.txt", "w"); fscanf(f, "%d %d", &npeaks, &np); - if (npeaks == 0) { + /*if (npeaks == 0) { fclose(g); fclose(f); remove("InitCondLS-temp.txt"); } - else { - fprintf(g, "%d %d\n", npeaks, np + ((npeaks > 0) ? 6 * nmod : 0)); + else {*/ + fprintf(g, "%d %d\n", npeaks, np + 6 * nmod ); peaks = (double*)malloc(sizeof(double) * npeaks); - printf("\nNumber of initial conditions: %d", np + ((npeaks > 0) ? 6 * nmod : 0)); + printf("\nNumber of initial conditions: %d", np + 6 * nmod); for (int i = 0; i < npeaks; i++) { fscanf(f, "%lg", &peaks[i]); fprintf(g, "%le", peaks[i]); @@ -1142,7 +1142,7 @@ int main(int argc, char* argv[]) { free(peaks); remove("InitCondLS.txt"); rename("InitCondLS-temp.txt", "InitCondLS.txt"); - } + //} } } else { From c0ad014f740922fad3371cd5477805e6c9845a54 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Mon, 23 Feb 2026 18:02:50 +0100 Subject: [PATCH 31/40] Add files via upload --- RTModel/include/LevMarFit.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RTModel/include/LevMarFit.h b/RTModel/include/LevMarFit.h index 1b3fa2c..45e4431 100644 --- a/RTModel/include/LevMarFit.h +++ b/RTModel/include/LevMarFit.h @@ -45,6 +45,8 @@ class LevMar { double tmaxmax; double maxmaxsum; + double y1maxmax; + double y2maxmax; double Tol; // void (LevMar::* PrintOut)(double*); From f657f68d191c576936730f35e4e98dce82957fa6 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Mon, 23 Feb 2026 19:25:02 +0100 Subject: [PATCH 32/40] Add files via upload --- RTModel/lib/LevMarFit.cpp | 8 ++++---- RTModel/lib/ModelSelector.cpp | 15 ++------------- 2 files changed, 6 insertions(+), 17 deletions(-) diff --git a/RTModel/lib/LevMarFit.cpp b/RTModel/lib/LevMarFit.cpp index 82812e5..b38e50f 100644 --- a/RTModel/lib/LevMarFit.cpp +++ b/RTModel/lib/LevMarFit.cpp @@ -321,8 +321,8 @@ void LevMar::ReadFiles(int argc, char* argv[]) { pr[1] = log(pr[1]); pr[4] = log(pr[4]); pr[5] = log(pr[5]); - /*current_path(exedir); - current_path("..");*/ + current_path(exedir); + current_path(".."); current_path("data"); VBM->LoadSunTable("SunEphemeris.txt"); current_path(eventname); @@ -1441,7 +1441,7 @@ double LevMar::ChiSquared(double* pr) { y1maxmax = y1max; y2maxmax = y2max; } - else if(fabs(maxsumn) > maxmaxsum); { + else if(fabs(maxsumn) > maxmaxsum) { maxmaxsum = - maxsumn; //calcolo del tempo || pr[6] è t0 in binary lens if (abs(double(pr[6] - t1)) > abs(double(pr[6] - t2))) { @@ -1470,7 +1470,7 @@ double LevMar::ChiSquared(double* pr) { y1maxmax = y1max; y2maxmax = y2max; } - else if (fabs(maxsumn) > maxmaxsum); { + else if (fabs(maxsumn) > maxmaxsum) { maxmaxsum = - maxsumn; if (abs(pr[6] - t1) > abs(pr[6] - t2)) { tmaxmax = t1; diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index 713e487..ac801c4 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -176,7 +176,6 @@ int main(int argc, char* argv[]) { supfac *= supfac; // Read curve to fit - printf("\n\nReading data\n"); current_path(eventname); @@ -291,7 +290,6 @@ int main(int argc, char* argv[]) { long end = ftell(f); - switch (modelcode[0]) { case 'P': sigmapr[1] = sigmapr[1] / pr[1]; @@ -337,13 +335,7 @@ int main(int argc, char* argv[]) { break; //da vedere per il triple lens - case 'T': - - - break; } - - if (c0 > 0) { if (nmod) { if (c0 < bumperlist->Amp) { @@ -609,7 +601,6 @@ int main(int argc, char* argv[]) { if (nmod > maxmodels) nmod = maxmodels; printf("\nModels to be saved = %d\n", nmod); - // Store best models passing the selection in directory "Models" current_path(".."); @@ -758,7 +749,6 @@ int main(int argc, char* argv[]) { } fclose(f); - scanbumper = bumperlist; for (il = 1; il <= nmod; il++) { for (int i = 0; i < nps; i++) { @@ -795,13 +785,12 @@ int main(int argc, char* argv[]) { xc = xc0; s0 = 0.5 * (sqrt(4 + xc * xc) - xc); - q2 = 0.00001; + q2 = 0.001; while (xc < 3 * sqrt(3 * q2) * s0 * s0 * s0) q2 *= 0.1; xc = xc0 - 3 * sqrt(3 * q2) * s0 * s0 * s0; s2 = 0.5 * (sqrt(4 + xc * xc) - xc); beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); - fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); @@ -948,7 +937,7 @@ int main(int argc, char* argv[]) { xc = xc0; s0 = 0.5 * (sqrt(4 + xc * xc) - xc); - q2 = 0.00001; + q2 = 0.001; while (xc < 3 * sqrt(3 * q2) * s0 * s0 * s0) q2 *= 0.1; xc = xc0 - 3 * sqrt(3 * q2) * s0 * s0 * s0; s2 = 0.5 * (sqrt(4 + xc * xc) - xc); From da5159b2b2afda772c10e545a3aa4ae36680d06d Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Mon, 23 Feb 2026 21:43:02 +0100 Subject: [PATCH 33/40] Add files via upload From 2650cb544b26d273fec97795fb38cd3ec4f425dc Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Mon, 23 Feb 2026 21:57:04 +0100 Subject: [PATCH 34/40] Add files via upload From aa52394cef18403f31de9bc4d8bd5fad011b2dc8 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Wed, 18 Mar 2026 15:07:00 +0100 Subject: [PATCH 35/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 91 +++++++++++++++++++++++++++++++++-- 1 file changed, 86 insertions(+), 5 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index ac801c4..9d7e131 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -176,6 +176,7 @@ int main(int argc, char* argv[]) { supfac *= supfac; // Read curve to fit + printf("\n\nReading data\n"); current_path(eventname); @@ -267,7 +268,7 @@ int main(int argc, char* argv[]) { for (int j = 0; j < nps + nlinpar * nfil; j++) { fscanf(f, "%le", &(pr[j])); } - + //per leggere i parametri fscanf(f, "%le", &(tmaxmax)); fscanf(f, "%le", &(y1maxmax)); @@ -290,6 +291,7 @@ int main(int argc, char* argv[]) { long end = ftell(f); + switch (modelcode[0]) { case 'P': sigmapr[1] = sigmapr[1] / pr[1]; @@ -335,7 +337,13 @@ int main(int argc, char* argv[]) { break; //da vedere per il triple lens + case 'T': + + + break; } + + if (c0 > 0) { if (nmod) { if (c0 < bumperlist->Amp) { @@ -601,6 +609,7 @@ int main(int argc, char* argv[]) { if (nmod > maxmodels) nmod = maxmodels; printf("\nModels to be saved = %d\n", nmod); + // Store best models passing the selection in directory "Models" current_path(".."); @@ -717,8 +726,8 @@ int main(int argc, char* argv[]) { } } } - - printf("\n- Adding initial conditions for planets"); + + printf("\n- Adding initial conditions for planets"); if (f = fopen("InitCondTS.txt", "r")) { int npeaks = 0; double* peaks; @@ -755,6 +764,77 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } + // Ci spostiamo nella directory "Models" + current_path(".."); + current_path("Models"); + + // Se 'filename' contiene un percorso (es. "Models/test.txt"), + // dobbiamo isolare solo "test.txt" perché siamo già dentro la cartella Models. + std::filesystem::path p(filename); + string solo_nome_file = p.filename().string(); + + // Tentiamo l'apertura usando il nome pulito + f = fopen(solo_nome_file.c_str(), "r"); + + if (f == NULL) continue; + // Un buffer abbastanza grande per contenere la prima riga + char riga[10000]; + + //Leggiamo solo la prima riga completa + if (fgets(riga, sizeof(riga), f) != NULL) { + + // cerchiamo i 5 double partendo dalla fine della riga poiché sappiamo che sono gli ultimi 5 valori della riga. + double t, y1, y2, m, c; // usiamo variabili temporanee + + // Cerchiamo di leggere i 5 valori. sscanf leggerà i primi 5 che trova se non specifichiamo altro, + + char* p = riga; + double val; + int count = 0; + vector tutti_i_numeri_della_riga; + + // Estraiamo tutti i numeri dalla riga per essere sicuri + char* endptr; + while (true) { + val = strtod(p, &endptr); + if (p == endptr) break; // Non ci sono più numeri + tutti_i_numeri_della_riga.push_back(val); + p = endptr; + } + + // Prendiamo gli ultimi 5 elementi trovati nella riga + size_t n = tutti_i_numeri_della_riga.size(); + if (n >= 5) { + c0 = tutti_i_numeri_della_riga[n - 1]; + maxmaxsum = tutti_i_numeri_della_riga[n - 2]; + y2maxmax = tutti_i_numeri_della_riga[n - 3]; + y1maxmax = tutti_i_numeri_della_riga[n - 4]; + tmaxmax = tutti_i_numeri_della_riga[n - 5]; + } + + } + //Assegnazione alla struttura bumperlist + bumperlist->tanomaly = tmaxmax; + bumperlist->y1anomaly = y1maxmax; + bumperlist->y2anomaly = y2maxmax; + bumperlist->maxsum = maxmaxsum; + bumperlist->Amp = c0; + + //printf("\nFile %s aperto e letto correttamente.", solo_nome_file.c_str()); + + /*else { + printf("\n! Errore: Formato dati non valido in %s", solo_nome_file.c_str()); + }*/ + fclose(f); + current_path(".."); + current_path("InitCond"); + //else { + // // Messaggio di debug per capire dove si trova il programma e cosa sta cercando + // printf("\n! Errore: Impossibile trovare %s in %s", + // solo_nome_file.c_str(), + // std::filesystem::current_path().string().c_str()); + //} + double u0 = pr[2]; double alpha = pr[3]; double s0, s = exp(pr[0]), s2; @@ -791,6 +871,7 @@ int main(int argc, char* argv[]) { s2 = 0.5 * (sqrt(4 + xc * xc) - xc); beta = beta0 + M_PI + asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); + fprintf(g, "%.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le %.10le\n", s, q, u0, alpha, rho, exp(pr[5]), pr[6], s2, q2, beta); beta = beta0 + M_PI - asin(fabs(2 * sqrt(q2 * (1 - s2 * s2)) / s2) / xc); @@ -937,7 +1018,7 @@ int main(int argc, char* argv[]) { xc = xc0; s0 = 0.5 * (sqrt(4 + xc * xc) - xc); - q2 = 0.001; + q2 = 0.00001; while (xc < 3 * sqrt(3 * q2) * s0 * s0 * s0) q2 *= 0.1; xc = xc0 - 3 * sqrt(3 * q2) * s0 * s0 * s0; s2 = 0.5 * (sqrt(4 + xc * xc) - xc); @@ -1450,4 +1531,4 @@ int main(int argc, char* argv[]) { free(Curv); return 0; -} +} \ No newline at end of file From a3ab90251a23af7d05368ff1ae0850a3107f8ddd Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Thu, 19 Mar 2026 11:28:31 +0100 Subject: [PATCH 36/40] Add files via upload --- RTModel/lib/ModelSelector.cpp | 128 ++++++++++++++++------------------ 1 file changed, 61 insertions(+), 67 deletions(-) diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index 9d7e131..78f49b0 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -235,7 +235,7 @@ int main(int argc, char* argv[]) { printf("\n- Model code: %s", modelcode); - + //da controllare pr = (double*)malloc(sizeof(double) * (nps + 20 + nlinpar * nfil)); //20 added to allow for higher order models in updates sigmapr = (double*)malloc(sizeof(double) * (nps + nlinpar * nfil)); Cov = (double*)malloc(sizeof(double) * (nps * nps)); @@ -280,12 +280,12 @@ int main(int argc, char* argv[]) { for (int j = 0; j < nps + nlinpar * nfil; j++) { fscanf(f, "%le", &(sigmapr[j])); - if (!(sigmapr[j] >= 0)) c0 = -1; + if (!(sigmapr[j] >= 0)) c0 = -1; } for (int i = 0; i < nps; i++) { for (int j = 0; j < nps; j++) { - fscanf(f, "%le", &(Cov[i + j * nps])); + fscanf(f, "%le", &(Cov[i + j * nps])); } } @@ -366,10 +366,10 @@ int main(int argc, char* argv[]) { scanbumper2->SetCovariance(Cov, dof / c0); strcpy(scanbumper2->modelcode, filename); scanbumper2->SetBuffer(f, start, end); - bumperlist->tanomaly = tmaxmax; - bumperlist->y1anomaly = y1maxmax; - bumperlist->y2anomaly = y2maxmax; - bumperlist->maxsum = maxmaxsum; + scanbumper2->tanomaly = tmaxmax; + scanbumper2->y1anomaly = y1maxmax; + scanbumper2->y2anomaly = y2maxmax; + scanbumper2->maxsum = maxmaxsum; scanbumper2->Amp = c0; scanbumper2->next = scanbumper->next; scanbumper->next = scanbumper2; @@ -764,77 +764,71 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } - // Ci spostiamo nella directory "Models" - current_path(".."); - current_path("Models"); + //// Ci spostiamo nella directory "Models" + //current_path(".."); + //current_path("Models"); - // Se 'filename' contiene un percorso (es. "Models/test.txt"), - // dobbiamo isolare solo "test.txt" perché siamo già dentro la cartella Models. - std::filesystem::path p(filename); - string solo_nome_file = p.filename().string(); + //// Se 'filename' contiene un percorso (es. "Models/test.txt"), + //// dobbiamo isolare solo "test.txt" perché siamo già dentro la cartella Models. + //std::filesystem::path p(filename); + //string solo_nome_file = p.filename().string(); - // Tentiamo l'apertura usando il nome pulito - f = fopen(solo_nome_file.c_str(), "r"); + //// Tentiamo l'apertura usando il nome pulito + //f = fopen(solo_nome_file.c_str(), "r"); - if (f == NULL) continue; - // Un buffer abbastanza grande per contenere la prima riga - char riga[10000]; + //if (f == NULL) continue; + //// Un buffer abbastanza grande per contenere la prima riga + //char riga[10000]; - //Leggiamo solo la prima riga completa - if (fgets(riga, sizeof(riga), f) != NULL) { + ////Leggiamo solo la prima riga completa + //if (fgets(riga, sizeof(riga), f) != NULL) { - // cerchiamo i 5 double partendo dalla fine della riga poiché sappiamo che sono gli ultimi 5 valori della riga. - double t, y1, y2, m, c; // usiamo variabili temporanee + // // cerchiamo i 5 double partendo dalla fine della riga poiché sappiamo che sono gli ultimi 5 valori della riga. + // double t, y1, y2, m, c; // usiamo variabili temporanee - // Cerchiamo di leggere i 5 valori. sscanf leggerà i primi 5 che trova se non specifichiamo altro, + // // Cerchiamo di leggere i 5 valori. sscanf leggerà i primi 5 che trova se non specifichiamo altro, - char* p = riga; - double val; - int count = 0; - vector tutti_i_numeri_della_riga; + // char* p = riga; + // double val; + // int count = 0; + // vector tutti_i_numeri_della_riga; - // Estraiamo tutti i numeri dalla riga per essere sicuri - char* endptr; - while (true) { - val = strtod(p, &endptr); - if (p == endptr) break; // Non ci sono più numeri - tutti_i_numeri_della_riga.push_back(val); - p = endptr; - } - - // Prendiamo gli ultimi 5 elementi trovati nella riga - size_t n = tutti_i_numeri_della_riga.size(); - if (n >= 5) { - c0 = tutti_i_numeri_della_riga[n - 1]; - maxmaxsum = tutti_i_numeri_della_riga[n - 2]; - y2maxmax = tutti_i_numeri_della_riga[n - 3]; - y1maxmax = tutti_i_numeri_della_riga[n - 4]; - tmaxmax = tutti_i_numeri_della_riga[n - 5]; - } + // // Estraiamo tutti i numeri dalla riga per essere sicuri + // char* endptr; + // while (true) { + // val = strtod(p, &endptr); + // if (p == endptr) break; // Non ci sono più numeri + // tutti_i_numeri_della_riga.push_back(val); + // p = endptr; + // } - } - //Assegnazione alla struttura bumperlist - bumperlist->tanomaly = tmaxmax; - bumperlist->y1anomaly = y1maxmax; - bumperlist->y2anomaly = y2maxmax; - bumperlist->maxsum = maxmaxsum; - bumperlist->Amp = c0; + // // Prendiamo gli ultimi 5 elementi trovati nella riga + // size_t n = tutti_i_numeri_della_riga.size(); + // if (n >= 5) { + // c0 = tutti_i_numeri_della_riga[n - 1]; + // maxmaxsum = tutti_i_numeri_della_riga[n - 2]; + // y2maxmax = tutti_i_numeri_della_riga[n - 3]; + // y1maxmax = tutti_i_numeri_della_riga[n - 4]; + // tmaxmax = tutti_i_numeri_della_riga[n - 5]; + // } - //printf("\nFile %s aperto e letto correttamente.", solo_nome_file.c_str()); - - /*else { - printf("\n! Errore: Formato dati non valido in %s", solo_nome_file.c_str()); - }*/ - fclose(f); - current_path(".."); - current_path("InitCond"); - //else { - // // Messaggio di debug per capire dove si trova il programma e cosa sta cercando - // printf("\n! Errore: Impossibile trovare %s in %s", - // solo_nome_file.c_str(), - // std::filesystem::current_path().string().c_str()); //} - + ////Assegnazione alla bumperlist + //bumperlist->tanomaly = tmaxmax; + //bumperlist->y1anomaly = y1maxmax; + //bumperlist->y2anomaly = y2maxmax; + //bumperlist->maxsum = maxmaxsum; + //bumperlist->Amp = c0; + + ////printf("\nFile %s aperto e letto correttamente.", solo_nome_file.c_str()); + + // /*else { + // printf("\n! Errore: Formato dati non valido in %s", solo_nome_file.c_str()); + // }*/ + //fclose(f); + //current_path(".."); + //current_path("InitCond"); + // double u0 = pr[2]; double alpha = pr[3]; double s0, s = exp(pr[0]), s2; From 7c8f15bfc910c20392e77c3e9a7657e3cf9ab0cb Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Thu, 26 Mar 2026 13:44:25 +0100 Subject: [PATCH 37/40] Add files via upload --- RTModel/lib/InitCond.cpp | 2 +- RTModel/lib/ModelSelector.cpp | 66 +---------------------------------- 2 files changed, 2 insertions(+), 66 deletions(-) diff --git a/RTModel/lib/InitCond.cpp b/RTModel/lib/InitCond.cpp index 560d17f..294c9a9 100644 --- a/RTModel/lib/InitCond.cpp +++ b/RTModel/lib/InitCond.cpp @@ -30,7 +30,7 @@ bool onlyorbital = false; // Only orbital motion models will be calculated. bool onlyupdate = false; // No model search, but only update of previously found best models. int usesatellite = 0; // Satellite to be used for initial conditions. Ground telescopes by default. char templatelibrary[256] = ""; // User-specified template library -char modelcategories[256] = "PSPXBSBOLSLXLOTXTS"; +char modelcategories[256] = "PSPXBSBOLSLXLOTSTX"; char astroini[256] = "0.0 0.0 0.125 1.0"; //double tau = 0.5; // Provisional!!! Exclude peaks shorter than tau diff --git a/RTModel/lib/ModelSelector.cpp b/RTModel/lib/ModelSelector.cpp index 78f49b0..77a2997 100644 --- a/RTModel/lib/ModelSelector.cpp +++ b/RTModel/lib/ModelSelector.cpp @@ -764,71 +764,7 @@ int main(int argc, char* argv[]) { pr[i] = scanbumper->p0[i]; } - //// Ci spostiamo nella directory "Models" - //current_path(".."); - //current_path("Models"); - - //// Se 'filename' contiene un percorso (es. "Models/test.txt"), - //// dobbiamo isolare solo "test.txt" perché siamo già dentro la cartella Models. - //std::filesystem::path p(filename); - //string solo_nome_file = p.filename().string(); - - //// Tentiamo l'apertura usando il nome pulito - //f = fopen(solo_nome_file.c_str(), "r"); - - //if (f == NULL) continue; - //// Un buffer abbastanza grande per contenere la prima riga - //char riga[10000]; - - ////Leggiamo solo la prima riga completa - //if (fgets(riga, sizeof(riga), f) != NULL) { - - // // cerchiamo i 5 double partendo dalla fine della riga poiché sappiamo che sono gli ultimi 5 valori della riga. - // double t, y1, y2, m, c; // usiamo variabili temporanee - - // // Cerchiamo di leggere i 5 valori. sscanf leggerà i primi 5 che trova se non specifichiamo altro, - - // char* p = riga; - // double val; - // int count = 0; - // vector tutti_i_numeri_della_riga; - - // // Estraiamo tutti i numeri dalla riga per essere sicuri - // char* endptr; - // while (true) { - // val = strtod(p, &endptr); - // if (p == endptr) break; // Non ci sono più numeri - // tutti_i_numeri_della_riga.push_back(val); - // p = endptr; - // } - - // // Prendiamo gli ultimi 5 elementi trovati nella riga - // size_t n = tutti_i_numeri_della_riga.size(); - // if (n >= 5) { - // c0 = tutti_i_numeri_della_riga[n - 1]; - // maxmaxsum = tutti_i_numeri_della_riga[n - 2]; - // y2maxmax = tutti_i_numeri_della_riga[n - 3]; - // y1maxmax = tutti_i_numeri_della_riga[n - 4]; - // tmaxmax = tutti_i_numeri_della_riga[n - 5]; - // } - - //} - ////Assegnazione alla bumperlist - //bumperlist->tanomaly = tmaxmax; - //bumperlist->y1anomaly = y1maxmax; - //bumperlist->y2anomaly = y2maxmax; - //bumperlist->maxsum = maxmaxsum; - //bumperlist->Amp = c0; - - ////printf("\nFile %s aperto e letto correttamente.", solo_nome_file.c_str()); - - // /*else { - // printf("\n! Errore: Formato dati non valido in %s", solo_nome_file.c_str()); - // }*/ - //fclose(f); - //current_path(".."); - //current_path("InitCond"); - // + double u0 = pr[2]; double alpha = pr[3]; double s0, s = exp(pr[0]), s2; From fadbae092b216bb515c36984cd7db213e54bc5f7 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Thu, 26 Mar 2026 13:44:55 +0100 Subject: [PATCH 38/40] Add files via upload From 51797aa8123970265a98fbe482932e5ae926849d Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Thu, 26 Mar 2026 14:43:47 +0100 Subject: [PATCH 39/40] Add files via upload From 2660edd1051d4277e2bf655a951cf82b8a05c396 Mon Sep 17 00:00:00 2001 From: Antonio Consiglio Date: Thu, 26 Mar 2026 14:49:37 +0100 Subject: [PATCH 40/40] Add files via upload