· 5 years ago · Oct 11, 2020, 08:54 PM
1/*
2 This file is part of the ULXLC model code.
3
4 The ULXLC model is free software: you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 any later version.
8
9 ULXLC is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 For a copy of the GNU General Public License see
14 <http://www.gnu.org/licenses/>.
15
16 Copyright 2016 Thomas Dauser, Remeis Observatory & ECAP
17
18 Modifcations 2020 by Norman Khan
19*/
20
21
22
23#include "ulxlcmod.h"
24#include <stdio.h>
25#include <stdlib.h>
26#include <time.h>
27#include <math.h>
28#include <sqlite3.h>
29#include <stdbool.h>
30
31/** global parameters, which can be used for several calls of the model */
32ulx_table* global_ulx_table=NULL;
33
34void check_par_model(inp_param* param){
35
36 // convert from degree to radian
37 double rad2deg = 180.0/PI;
38
39 if (param->beta<0){
40 printf(" *** warning, value of beta=%e < 0 not allowed. Resetting to beta=0.0 \n",param->beta);
41 param->beta=0.0;
42 }
43 if (param->beta>0.999){
44 printf(" *** warning, value of beta=%e > 0.999 not allowed. Resetting to beta=0.999 \n",param->beta);
45 param->beta=0.999;
46 }
47
48 if (param->theta<0){
49 printf(" *** warning, value of theta=%e < 2 deg not allowed. Resetting to theta=0.0 \n",param->theta*rad2deg);
50 param->theta=2.0/rad2deg;
51 }
52 if (param->theta>PI/4.0){
53 printf(" *** warning, value of theta=%e > 45.0 deg not allowed. Resetting to theta=45.0 deg \n",param->theta*rad2deg);
54 param->theta=PI/4.0;
55 }
56
57 if (param->incl<0){
58 printf(" *** warning, value of incl=%e < 0 deg not allowed. Resetting to incl=0.0 deg \n",param->incl*rad2deg);
59 param->incl=0.0;
60 }
61 if (param->incl>PI/2.0){
62 printf(" *** warning, value of incl=%e > 90 deg not allowed. Resetting to incl=90 deg \n",param->incl*rad2deg);
63 param->incl=PI/2.0;
64 }
65 if (param->dopulse>1 || param->dopulse<0){
66 printf(" *** warning, value of dopulse = %i, but can only be 0 (light curve) or 1 (pulse) \n",param->dopulse);
67 param->dopulse=0;
68 }
69
70 return;
71}
72
73void init_par_model(inp_param* param, const double* inp_par, const int n_parameter, int* status){
74
75 assert(n_parameter == NUM_PARAM);
76
77 // convert from degree to radian
78 double deg2rad = PI/180.0;
79
80 param->period = inp_par[0];
81 param->phase = inp_par[1];
82 param->theta = inp_par[2]*deg2rad;
83 param->incl = inp_par[3]*deg2rad;
84 param->dincl = inp_par[4]*deg2rad;
85 param->beta = inp_par[5];
86 param->dopulse = (int) inp_par[6];
87 check_par_model(param);
88}
89
90// calculate the flux array for the current value of theta
91static double* interp_betTheta(ulx_table* tab,double theta, double beta, int* status){
92
93 assert(tab!=NULL);
94
95 int ind_th = binary_search(theta, tab->theta, tab->ntheta);
96 if (ind_th == -1){
97 printf(" *** error: values of theta=%.2f not tabulated",theta*180/3.1415);
98 ULXLC_ERROR(" failed to look up the value of theta in the loaded table ",status);
99 return NULL;
100 }
101
102 int ind_bet = binary_search(beta, tab->beta, tab->nbeta);
103 if (ind_bet == -1){
104 printf(" *** error: values of beta=%e not tabulated",beta);
105 ULXLC_ERROR(" failed to look up the value of beta in the loaded table ",status);
106 return NULL;
107 }
108
109 // allocate memory for the flux array
110 double* inp_flux = (double*) malloc (tab->nang*sizeof(double));
111 CHECK_MALLOC_RET_STATUS(inp_flux,status,NULL);
112
113 // 2d-interpolate the flux now
114 double fac_th = (theta-tab->theta[ind_th]) / (tab->theta[ind_th+1]-tab->theta[ind_th]);
115 double fac_bet = (beta-tab->beta[ind_bet]) / (tab->beta[ind_bet+1]-tab->beta[ind_bet]);
116 int ii;
117 for (ii=0;ii<tab->nang;ii++){
118 inp_flux[ii] = (1-fac_bet)*(1-fac_th) * tab->flux[ind_bet][ind_th][ii]
119 + (fac_bet)*(1-fac_th) * tab->flux[ind_bet+1][ind_th][ii]
120 + (1-fac_bet)*(fac_th) * tab->flux[ind_bet][ind_th+1][ii]
121 + (fac_bet)*(fac_th) * tab->flux[ind_bet+1][ind_th+1][ii];
122 }
123
124 return inp_flux;
125}
126
127static double calc_orbit_ang(double phase, inp_param param){
128 // need the sine distributed from [0,1]
129 double ang = 0.5*(sin(phase*2*PI+PI/2)+1)*(param.dincl*2) + param.incl - param.dincl;
130 return fabs(ang); // return the absolute value
131}
132
133static double get_flux_from_angle(double ang, double* flux, double* ang_arr, int nang, int* status){
134
135 double ang_min = 0.0;
136 double ang_max = PI/2;
137
138 // allow angles between 90 to 180 degrees
139 if ( ang > PI/2 && ang < PI){
140 ang = PI - ang;
141 }
142
143 // check at the boundaries
144 if (ang < ang_arr[0] && ang>=ang_min ){
145 return flux[0];
146 } else if (ang > ang_arr[nang-1] && ang <= ang_max){
147 return flux[nang-1];
148 }
149
150 // get the corresponding index
151 int ind = binary_search(ang,ang_arr,nang);
152
153 // check if something went wrong (should not happen)
154 if (ind == -1){
155 printf(" *** error: angle %e [deg] is not tabulated (shoud be in [%.2e,%.2e]\n",ang*180/PI,
156 ang_arr[0]*180/PI,ang_arr[nang-1]*180/PI);
157 ULXLC_ERROR("failed to get the flux for the desired phase",status);
158 return -1.0;
159 }
160
161 // interpolate now the flux
162 double fac = (ang-ang_arr[ind]) / (ang_arr[ind+1]-ang_arr[ind]);
163 return (1-fac)*flux[ind]+fac*flux[ind+1];
164}
165
166// input: t[nt+1], photar[nt]
167
168static void create_lc(const double* t0,const int nt,double* photar,
169 double* flux, ulx_table* tab,inp_param param,int* status){
170
171 int ii;
172 for (ii=0; ii<nt; ii++){ //nt as photar[nt]
173 // if we want only a single pulse
174
175 double phase=0.0;
176
177 if (param.dopulse==1){
178 // take the phase in the middle again
179 double delta_t_phase_lo = (t0[ii] - t0[0])/ param.period - param.phase ;
180 double delta_t_phase_hi = (t0[ii+1] - t0[0])/ param.period - param.phase ;
181 phase=0.5*(delta_t_phase_lo+delta_t_phase_hi);
182 if ( ( phase < 0 ) || (phase >= 1 ) ){
183 continue;
184 }
185
186 } else {
187
188 // see which part in the period we are in (then time is between 0 and param.period)
189 phase = (0.5*(t0[ii]+t0[ii+1])/param.period + param.phase) ;
190 phase = phase - (double) ((int) phase); // just get everything behind the colon
191
192 }
193
194
195
196 // now get the corresponding angle
197 double ang = calc_orbit_ang(phase,param);
198 photar[ii] = get_flux_from_angle(ang,flux,tab->ang,tab->nang,status);
199 // printf("t: %d \n", ii);
200 // printf("flux: %f \n", flux);
201 // printf("tab->ang %lf \n", tab->ang);
202 // printf("tab->nang %lf \n", tab->nang);
203 // printf("Photar[t]: %f \n", photar[ii]);
204 // printf("status: %d \n", *status);
205 CHECK_STATUS_VOID(*status);
206 }
207
208 return;
209}
210
211
212// basic calculation of the model
213// input: t[nt+1], photar[nt]
214// output: photar
215static void ulxlc_base(const double* t, const int nt, double* photar, inp_param param_struct,int* status){
216 // load the GLOBAL TABLE
217 if (global_ulx_table==NULL){
218 ulxlc_load_table(ULXTABLE_FILENAME, &global_ulx_table, status);
219 CHECK_ULXLC_ERROR("reading the table failed",status);
220 CHECK_STATUS_VOID(*status);
221 }
222 assert(global_ulx_table!=NULL);
223
224 // interpolate the flux for a given opening_angle/theta
225 double* flux = interp_betTheta(global_ulx_table, param_struct.theta,param_struct.beta,status);
226 CHECK_ULXLC_ERROR("reading of theta values from table failed",status);
227 CHECK_STATUS_VOID(*status);
228 create_lc(t,nt,photar, flux,global_ulx_table,param_struct,status);
229 free(flux);
230 CHECK_STATUS_VOID(*status);
231
232// free_ulxTable(&global_ulx_table);
233}
234
235void free_ulxTable(void){
236 ulx_table* tab = global_ulx_table;
237 if (tab!=NULL){
238 free(tab->theta);
239 free(tab->ang);
240 free(tab->beta);
241
242 if(tab->flux!=NULL){
243 int ii,jj;
244 for (jj=0; jj<tab->nbeta; jj++){
245
246 if (tab->flux[jj]!=NULL){
247 for(ii=0;ii<tab->ntheta;ii++){
248 free(tab->flux[jj][ii]);
249 }
250 free(tab->flux[jj]);
251 }
252
253 }
254 free(tab->flux);
255 }
256 free(tab);
257 }
258 // set it to NULL to know the table is not there anymore
259 tab = NULL;
260}
261
262int ulxlc_model(const double* t, const int nt, double* photar, const double* parameter, const int n_parameter){
263
264 int status = EXIT_SUCCESS;
265 // fill in parameters
266 inp_param param_struct;
267
268 init_par_model(¶m_struct,parameter,n_parameter,&status);
269 CHECK_STATUS_RET(status,status);
270
271 // call the function which calculates the light curve
272 ulxlc_base(t, nt, photar, param_struct,&status);
273 CHECK_STATUS_RET(status,status);
274
275 return status;
276}
277
278
279/*
280##########################
281Added functions start here
282##########################
283*/
284
285// Auxilary (Things that any reasonable programming language should do)
286double max(double* arr, int len){
287 // Get the maximum value of an array
288 // arr : array
289 // len : length of array
290 int i;
291 double t;
292 t = arr[0];
293 for(i=1; i<len ;i++)
294 {
295 if(arr[i]>t)
296 t=arr[i];
297 }
298 return(t);
299}
300
301double min(double* arr, int len){
302 // Get the minimum value of an array
303 // x : array
304 // len : length of array
305 int i;
306 double t;
307 t = arr[0];
308 for(i=1; i<len ;i++)
309 {
310 if(arr[i]<t)
311 t=arr[i];
312 }
313 return(t);
314}
315
316
317int classify_curve(double lc_ulx_lim, double lc_max_flux, double lc_min_flux){
318 // 0 : dead
319 // 1 : transient
320 // 2 : alive
321 int classification;
322
323 if (lc_min_flux > lc_ulx_lim){
324 classification = 2;
325 }
326 else if (lc_max_flux < lc_ulx_lim){
327 classification = 0;
328 }
329 else if (lc_ulx_lim < lc_max_flux && lc_ulx_lim > lc_min_flux){
330 classification = 1;
331 }
332 return classification;
333}
334
335int get_first_transient_cycle_index(int is_ulx[]){
336 int diff;
337 for(int i=1; i<8; i++){
338 diff = is_ulx[i] - is_ulx[i-1];
339 if (diff != 0){
340 return i;
341 }
342 }
343 return -1;
344}
345
346
347
348// some sort of sqlite thingy idk
349static int callback(void *data, int argc, char **argv, char **azColName) {
350 int i;
351
352 fprintf(stderr, "%s: \n", (const char*)data);
353
354 for(i = 0; i<argc; i++) {
355 printf("%d \n", i);
356 printf("%s = %s\n", azColName[i], argv[i] ? argv[i] : "NULL");
357 }
358 printf("\n");
359 return 0;
360}
361
362
363int main(int argc, char *argv[]){
364 // argv[1] = run_id
365 // argv[2] = size (system_N)
366 // argv[3] = erass_system_period_cutoff
367 // argv[4] = lmxrb_duty_cycle
368
369 if(argc!=5){
370 printf("Incorrect number of parameters passed (%d / %d) \n", argc-1, 4);
371 exit(1);
372 }
373
374 double fmod(double x, double y); // Modulo Operator
375 srand(time(NULL)); // Initialise random seed
376
377 char* run_id = argv[1]; // run identifier
378
379 // ULXLC parameters
380 double params[7];
381 params[0] = 50.0; // period
382 params[1] = 0.0; // phase
383 params[2] = NAN; // theta
384 params[3] = NAN; // incl
385 params[4] = NAN; // dincl
386 params[5] = 0.3; // beta
387 params[6] = 0; // dopulse
388
389 // Setup Curve
390 int lc_nt = 5000; // Length of time series
391 float lc_timestep = 0.01; // Timestep in seconds
392 float lc_duration = lc_timestep*lc_nt; // Lightcurve duration
393 double lc_t[lc_nt]; // Time array
394 double lc_flux[lc_nt]; // Photon detection Array (Flux)
395
396 // lightcurve properties
397 double lc_period = params[0];
398 double lc_zero_incl_max_flux;
399 double lc_max_flux;
400 double lc_min_flux;
401 double lc_flux_scaling_constant;
402 double lc_P_wind_time_scaling_constant;
403 double lc_P_sup_time_scaling_constant;
404 double lc_ulx_lim=0;
405 int lc_classification; // 0 = Dead, 1 = Transient, 2 = Dead
406
407 // System parameters
408 int system_N = atoi(argv[2]);
409 int system_row_id[system_N];
410 double system_theta[system_N];
411 int system_inclination[system_N];
412 int system_dincl[system_N];
413 double system_Lx[system_N]; // In units of 10^39 erg s^-1
414 double system_P_wind_days[system_N];
415 double system_P_sup_days[system_N];
416 int system_lmxrb[system_N]; // Is the system defined as a LMXRB that could go into outburst?
417
418 // Not Used
419 int system_mttype = 0; // Mass transfer type: 0 = None, 1 = Nuclear, 2 = Thermal, 3 = WD mass transfer
420
421 // eRASS Curve Sampling
422 int erass_system_period_cutoff = atoi(argv[3]); // Treat sources that are over a wind period length as persistent sources
423 int erass_sample_interval = 30*6; // Sampling interval in days
424 int erass_sample_iterations = 10000; // Number of Monte Carlo cycle repeats
425
426 float erass_lmxrb_duty_cycle = atof(argv[4]); // Duty Cycle for LMXRBs used in eRASS simulation
427
428 int erass_P_wind_start_time_days[erass_sample_iterations]; // Used to hold time of first observation
429 int erass_P_sup_start_time_days[erass_sample_iterations]; // Used to hold time of first observation
430
431 int erass_P_wind_sample_time_days[erass_sample_iterations][8]; // Sample times in days
432 int erass_P_sup_sample_time_days[erass_sample_iterations][8]; // Sample times in days
433
434 double erass_P_wind_lc_sample_times[erass_sample_iterations][8]; // Sample times on the curve
435 double erass_P_sup_lc_sample_times[erass_sample_iterations][8]; // Sample times on the curve
436
437 int erass_P_wind_lc_sample_indexs[erass_sample_iterations][8]; // Sample time index on the curve
438 int erass_P_sup_lc_sample_indexs[erass_sample_iterations][8]; // Sample time index on the curve
439
440 double erass_P_wind_lc_sample_fluxes[erass_sample_iterations][8]; // Sample fluxes on the curve
441 double erass_P_sup_lc_sample_fluxes[erass_sample_iterations][8]; // Sample fluxes on the curve
442
443 double erass_P_wind_sample_fluxes[erass_sample_iterations][8]; // Sample fluxes in x10^39 erg s^-1
444 double erass_P_sup_sample_fluxes[erass_sample_iterations][8]; // Sample fluxes in x10^39 erg s^-1
445
446 int erass_P_wind_is_ulx[erass_sample_iterations][8]; // Is the sampled time observed as a ULX?
447 int erass_P_sup_is_ulx[erass_sample_iterations][8]; // Is the sampled time observed as a ULX?
448
449 double erass_1_ulx_prob; // Probability of source being a ulx on first erass cycle
450
451 double erass_P_wind_persistent_prob; // probability of the source being persistent for the entire length of eRASS
452 double erass_P_sup_persistent_prob; // probability of the source being persistent for the entire length of eRASS
453
454 double erass_P_wind_transient_prob[8]; // probability of the source being transient for the entire length of eRASS
455 double erass_P_sup_transient_prob[8]; // probability of the source being transient for the entire length of eRASS
456
457
458 // SQLite
459 char filename_database[9];
460 sprintf(filename_database, "ulxlc.db");
461 sqlite3 *db;
462 char *zErrMsg = 0;
463 char *sql;
464 int rc;
465 const char* data = "Callback function called";
466 sqlite3_stmt *res;
467
468 rc = sqlite3_open(filename_database, &db);
469
470 if(rc) {
471 fprintf(stderr, "Can't open database: %s\n", sqlite3_errmsg(db));
472 return(0);
473 } else {
474 fprintf(stderr, "Opened database successfully\n");
475 }
476
477 // Create CLASSIFICATIONS SQL table
478 sql = "CREATE TABLE IF NOT EXISTS CLASSIFICATIONS(" \
479 "system_row_id INT," \
480 "system_theta REAL," \
481 "system_dincl INT," \
482 "system_inclination INT," \
483 "lc_min_flux REAL," \
484 "lc_max_flux REAL," \
485 "lc_ulx_lim REAL," \
486 "lc_classification INT," \
487 "run_id CHAR);";
488
489 rc = sqlite3_exec(db, sql, callback, 0, &zErrMsg);
490 if( rc != SQLITE_OK ){
491 fprintf(stderr, "SQL error: %s\n", zErrMsg);
492 sqlite3_free(zErrMsg);
493 } else {
494 fprintf(stdout, "CLASSIFICATIONS Table created successfully\n");
495 }
496
497 // Create Transient results SQL table
498 sql = "CREATE TABLE IF NOT EXISTS TRANSIENT(" \
499 "system_row_id INT," \
500 "system_dincl INT," \
501 "system_inclination INT," \
502 "erass_1_ulx_prob REAL," \
503 "erass_2_P_wind_transient_prob REAL," \
504 "erass_3_P_wind_transient_prob REAL," \
505 "erass_4_P_wind_transient_prob REAL," \
506 "erass_5_P_wind_transient_prob REAL," \
507 "erass_6_P_wind_transient_prob REAL," \
508 "erass_7_P_wind_transient_prob REAL," \
509 "erass_8_P_wind_transient_prob REAL," \
510 "erass_2_P_sup_transient_prob REAL," \
511 "erass_3_P_sup_transient_prob REAL," \
512 "erass_4_P_sup_transient_prob REAL," \
513 "erass_5_P_sup_transient_prob REAL," \
514 "erass_6_P_sup_transient_prob REAL," \
515 "erass_7_P_sup_transient_prob REAL," \
516 "erass_8_P_sup_transient_prob REAL," \
517 "erass_P_wind_persistent_prob REAL," \
518 "erass_P_sup_persistent_prob REAL," \
519 "run_id CHAR);";
520
521
522
523
524 rc = sqlite3_exec(db, sql, callback, 0, &zErrMsg);
525 if( rc != SQLITE_OK ){
526 fprintf(stderr, "SQL error: %s\n", zErrMsg);
527 sqlite3_free(zErrMsg);
528 } else {
529 fprintf(stdout, "TRANSIENT Table created successfully\n");
530 }
531
532
533 /* Retrieve sim rows from sql table */
534 sql = sqlite3_mprintf("SELECT row_id, theta_half_deg, Lx1, P_wind_days, P_sup_days, inclination, dincl, lmxrb from ERASS_MC_SAMPLED_SYSTEMS WHERE run_id=%Q", run_id);
535 // sql = sqlite3_mprintf("SELECT * from ERASS_MC_SAMPLED_SYSTEMS WHERE run_id='%Q'", run_id);
536 rc = sqlite3_open(filename_database, &db);
537 rc = sqlite3_prepare_v2(db, sql, -1, &res, 0);
538 int r = 0;
539 while(sqlite3_step(res) == SQLITE_ROW){
540 system_row_id[r] = sqlite3_column_int(res, 0);
541 system_theta[r] = sqlite3_column_double(res, 1);
542 system_Lx[r] = sqlite3_column_double(res, 2);
543 system_P_wind_days[r] = sqlite3_column_double(res, 3);
544 system_P_sup_days[r] = sqlite3_column_double(res, 4);
545 system_inclination[r] = sqlite3_column_int(res, 5);
546 system_dincl[r] = sqlite3_column_int(res, 6);
547 system_lmxrb[r] = sqlite3_column_int(res, 7);
548 r++;
549 }
550 sqlite3_finalize(res);
551 sqlite3_close(db);
552
553 // Populate time array
554 for (int i=0; i<lc_nt; i++){
555 lc_t[i] = lc_timestep*i;
556 }
557
558
559
560 for (int N=0;N<system_N;N++){
561 params[2] = system_theta[N]; // theta
562 params[3] = 0; // incl
563 params[4] = system_dincl[N]; // dincl
564 if (system_theta[N] > 45){
565 continue;
566 }
567 // Run the lightcurve model.
568 ulxlc_model(lc_t, lc_nt, lc_flux, params, 7);
569
570 // 0 Inclination is treated as maximum possible Lx (Looking straight down the windcone)
571 if (params[3]==0){
572 // Calculate Flux normalisation
573 lc_zero_incl_max_flux = max(lc_flux, lc_nt);
574 lc_flux_scaling_constant = system_Lx[N] / lc_zero_incl_max_flux; // Curve Normalisation constant
575 lc_ulx_lim = 1 / lc_flux_scaling_constant; //Units of 10^39 erg s^-1
576 }
577
578 params[3] = system_inclination[N]; // incl
579
580 // Run the lightcurve model.
581 ulxlc_model(lc_t, lc_nt, lc_flux, params, 7);
582
583 // Classify the lightcurve as alive/dead/transient
584 lc_max_flux = max(lc_flux, lc_nt);
585 lc_min_flux = min(lc_flux, lc_nt);
586 lc_classification = classify_curve(lc_ulx_lim, lc_max_flux, lc_min_flux);
587
588 // TODO
589 int erass_sample = 1; // Sample erass yes/no
590
591 if (system_lmxrb[N] == 1){
592 // printf("System could be in outburst \n");
593 double r = ((double) rand() / (RAND_MAX));
594 // printf("r = %f | duty cycle = %f\n", r, erass_lmxrb_duty_cycle);
595 if (r > erass_lmxrb_duty_cycle){
596 int erass_sample = 0;
597 // printf("NOT SAMPLING zn");
598 }
599
600 }
601 // Lightcurve is transient & P_wind or P_sup is below limit
602 if ((lc_classification == 1) & (erass_sample==1) & (system_P_wind_days[N] < erass_system_period_cutoff || system_P_sup_days[N] < erass_system_period_cutoff) ){
603
604 int erass_1_N_ulx = 0; // Number of times the system was observed as a ULX on the first cycle.
605
606 int erass_P_wind_N_transient[8] = {0,0,0,0,0,0,0,0}; // Probability of lc being identified as transient at cycle i
607 int erass_P_sup_N_transient[8] = {0,0,0,0,0,0,0,0}; // Probability of lc being identified as transient at cycle i
608
609 int erass_P_wind_N_persistent = 0; // Times lc was being identified as persistent over eRASS
610 int erass_P_sup_N_persistent = 0; // Times lc was being identified as persistent over eRASS
611
612
613 for (int j=0; j<erass_sample_iterations; j++){ // MC lightcurve sampling repeats
614
615 // Calculate scaling constants
616 lc_P_wind_time_scaling_constant = lc_period / system_P_wind_days[N]; // Time scaling constant
617 lc_P_sup_time_scaling_constant = lc_period / system_P_sup_days[N]; // Time scaling constant
618
619 // Sampling in days
620 erass_P_wind_start_time_days[j] = (rand() % ((int) system_P_wind_days[N] - 0 + 1)) + 0; // Random integer between 0 and P_wind_days
621 erass_P_sup_start_time_days[j] = (rand() % ((int) system_P_sup_days[N] - 0 + 1)) + 0; // Random integer between 0 and system_P_sup_days
622
623 // printf("eRASS Sampling: j = %d / %d \n", j, erass_sample_iterations);
624 // printf("------------------------------- \n");
625 // printf("erass_system_period_cutoff: %d days \n", erass_system_period_cutoff);
626 // printf("erass_sample_interval: %d days \n", erass_sample_interval);
627 // printf("erass_P_wind_start_time_days: %d days \n", erass_P_wind_start_time_days[j]);
628 // printf("erass_P_sup_start_time_days: %d days \n", erass_P_sup_start_time_days[j]);
629 // printf("system_P_wind_days %f \n", system_P_wind_days[N]);
630 // printf("system_P_sup_days %f \n", system_P_sup_days[N]);
631 // printf("lc_P_wind_time_scaling_constant %f \n", lc_P_wind_time_scaling_constant);
632 // printf("lc_P_sup_time_scaling_constant %f \n", lc_P_sup_time_scaling_constant);
633 // printf("---------------------------------------------------------------------------------------------------------------- \n");
634
635
636 for (int i=0; i<8; i++){
637 erass_P_wind_sample_time_days[j][i] = erass_P_wind_start_time_days[j] + i * erass_sample_interval;
638 erass_P_sup_sample_time_days[j][i] = erass_P_sup_start_time_days[j] + i * erass_sample_interval;
639
640 erass_P_wind_lc_sample_times[j][i] = fmod(lc_P_wind_time_scaling_constant *erass_P_wind_sample_time_days[j][i], lc_period); // Sampling points scaled for curve
641 erass_P_sup_lc_sample_times[j][i] = fmod(lc_P_sup_time_scaling_constant * erass_P_sup_sample_time_days[j][i], lc_period); // Sampling points scaled for curve
642
643 erass_P_wind_lc_sample_indexs[j][i] = binary_search(erass_P_wind_lc_sample_times[j][i], lc_t, lc_nt); // Get corresponding indexs
644 erass_P_sup_lc_sample_indexs[j][i] = binary_search(erass_P_sup_lc_sample_times[j][i], lc_t, lc_nt); // Get corresponding indexs
645
646 erass_P_wind_lc_sample_fluxes[j][i] = lc_flux[erass_P_wind_lc_sample_indexs[j][i]]; // Find corresponding fluxes on curve
647 erass_P_sup_lc_sample_fluxes[j][i] = lc_flux[erass_P_sup_lc_sample_indexs[j][i]]; // Find corresponding fluxes on curve
648
649 // erass_P_wind_sample_fluxes[j][i] = lc_flux_scaling_constant * erass_P_wind_lc_sample_fluxes[j][i]; // Fluxes in 1x10^39 erg s^-1 units
650 // erass_P_sup_sample_fluxes[j][i] = lc_flux_scaling_constant * erass_P_sup_lc_sample_fluxes[j][i]; // Fluxes in 1x10^39 erg s^-1 units
651
652 erass_P_wind_is_ulx[j][i] = erass_P_wind_lc_sample_fluxes[j][i] > lc_ulx_lim; // True/False is sampled flux ulx?
653 erass_P_sup_is_ulx[j][i] = erass_P_sup_lc_sample_fluxes[j][i] > lc_ulx_lim; // True/False is sampled flux ulx?
654
655 // printf("i=%d \n", i);
656 // printf("erass_P_wind_start_time_days: %d days | erass_P_sup_start_time_days: %d days \n", erass_P_wind_start_time_days[j], erass_P_sup_start_time_days[j]);
657 // printf("erass_P_wind_sample_time_days[i]: %d days | erass_P_sup_sample_time_days[i]: %d days \n", erass_P_wind_sample_time_days[j][i], erass_P_sup_sample_time_days[j][i]);
658 // printf("erass_P_wind_lc_sample_times[i]: %f units | erass_P_sup_lc_sample_times[i]: %f units \n", erass_P_wind_lc_sample_times[j][i], erass_P_sup_lc_sample_times[j][i]);
659 // printf("erass_P_wind_lc_sample_indexs[i]: %d | erass_P_sup_lc_sample_indexs[i]: %d \n", erass_P_wind_lc_sample_indexs[j][i], erass_P_sup_lc_sample_indexs[j][i]);
660 // printf("erass_P_wind_lc_sample_fluxes[i]: %f | erass_P_sup_lc_sample_fluxes[i]: %f\n", erass_P_wind_lc_sample_fluxes[j][i], erass_P_sup_lc_sample_fluxes[j][i]);
661 // printf("erass_P_wind_sample_fluxes[i]: %f x10^39 erg s^-1 | erass_P_sup_sample_fluxes[i]: %f x10^39 erg s^-1\n", erass_P_wind_sample_fluxes[j][i], erass_P_sup_sample_fluxes[j][i]);
662 // printf("erass_P_wind_is_ulx[i] %d | erass_P_sup_is_ulx[i] %d \n", erass_P_wind_is_ulx[j][i], erass_P_sup_is_ulx[j][i]);
663 // printf("---------------------------------------------------------------------------------------------------------------- \n");
664 }
665 // printf("================================================================================================================ \n");
666
667
668
669 erass_1_N_ulx = erass_1_N_ulx + erass_P_wind_is_ulx[j][0]; // This probability is equal for both P_wind and P_orb as uniform sampling in the first cycle is irrelevant of length.
670
671 int ftci_P_wind;
672 int ftci_P_sup;
673
674 ftci_P_wind = get_first_transient_cycle_index(erass_P_wind_is_ulx[j]); // Get the first cycle for which the source was observed as transent
675 ftci_P_sup = get_first_transient_cycle_index(erass_P_sup_is_ulx[j]); // Get the first cycle for which the source was observed as transent
676
677 // Calculate final probabilities
678 if (system_P_wind_days[N] < erass_system_period_cutoff){
679 if (ftci_P_wind==-1){
680 erass_P_wind_N_persistent++;
681 }
682 else{
683 erass_P_wind_N_transient[ftci_P_wind]++;
684 }
685 }
686 else{
687 erass_P_wind_N_persistent = erass_sample_iterations; // let every observation be a persistent one if system_P_wind_days > erass_system_period_cutoff
688 }
689
690 if (system_P_sup_days[N] < erass_system_period_cutoff){
691 if (ftci_P_sup==-1){
692 erass_P_sup_N_persistent++;
693 }
694 else{
695 erass_P_sup_N_transient[ftci_P_sup]++;
696 }
697 // printf("%d \n", erass_P_wind_N_transient[1]);
698 }
699 else{
700 erass_P_sup_N_persistent = erass_sample_iterations; // let every observation be a persistent one if system_P_sup_days > erass_system_period_cutoff
701 }
702
703 }
704
705 erass_1_ulx_prob = (double) erass_1_N_ulx / erass_sample_iterations;
706 erass_P_wind_persistent_prob = (double) erass_P_wind_N_persistent / erass_sample_iterations;
707 erass_P_sup_persistent_prob = (double) erass_P_sup_N_persistent / erass_sample_iterations;
708
709 for (int i=0; i<8;i++){
710 erass_P_wind_transient_prob[i] = (double) erass_P_wind_N_transient[i] / erass_sample_iterations;
711 erass_P_sup_transient_prob[i] = (double) erass_P_sup_N_transient[i] / erass_sample_iterations;
712 }
713
714
715
716 // Insert TRANSIENT into sqlite table
717 char* sql2 = sqlite3_mprintf("INSERT INTO TRANSIENT VALUES (" \
718 "%d, %d, %d, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %Q);",
719 system_row_id[N],
720 system_dincl[N],
721 system_inclination[N],
722 erass_1_ulx_prob,
723 erass_P_wind_transient_prob[1],
724 erass_P_wind_transient_prob[2],
725 erass_P_wind_transient_prob[3],
726 erass_P_wind_transient_prob[4],
727 erass_P_wind_transient_prob[5],
728 erass_P_wind_transient_prob[6],
729 erass_P_wind_transient_prob[7],
730 erass_P_sup_transient_prob[1],
731 erass_P_sup_transient_prob[2],
732 erass_P_sup_transient_prob[3],
733 erass_P_sup_transient_prob[4],
734 erass_P_sup_transient_prob[5],
735 erass_P_sup_transient_prob[6],
736 erass_P_sup_transient_prob[7],
737 erass_P_wind_persistent_prob,
738 erass_P_sup_persistent_prob,
739 run_id);
740
741 rc = sqlite3_open(filename_database, &db);
742 rc = sqlite3_exec(db, sql2, callback, 0, &zErrMsg);
743
744 if( rc != SQLITE_OK ){
745 fprintf(stderr, "SQL error: %s\n", zErrMsg);
746 sqlite3_free(zErrMsg);
747 } else {
748 // fprintf(stdout, "TRANSIENT Records created successfully\n");
749 }
750 sqlite3_close(db);
751 }
752
753 // printf("\n");
754 // printf("System Information: (ID %d) (iter = %d) \n", system_row_id, N);
755 // printf("--------------------------------------- \n");
756 // printf("system_theta: %f deg\n", system_theta);
757 // printf("system_dincl: %d deg\n", system_dincl);
758 // printf("system_inclination: %d deg\n", system_inclination);
759 // printf("system_Lx: %f x 10^39 erg s^-1 \n", system_Lx);
760 // printf("system_P_wind_days: %f days\n", system_P_wind_days);
761 // printf("system_P_sup_days: %f days\n", system_P_sup_days);
762 // printf("system_mttype: %d \n", system_mttype[0]);
763 // printf("system_lmxrb: %f \n", system_lmxrb[0]);
764 //
765 // printf("\n");
766 //
767 // printf("Lightcurve Information: \n");
768 // printf("----------------------- \n");
769 // printf("lc_nt: %d \n", lc_nt);
770 // printf("lc_timestep: %f seconds \n", lc_timestep);
771 // printf("lc_duration: %f seconds \n", lc_duration);
772 // printf("lc_period: %f seconds \n", lc_period);
773 // printf("lc_zero_incl_max_flux: %f units\n", lc_zero_incl_max_flux);
774 // printf("lc_max_flux: %f units\n", lc_max_flux);
775 // printf("lc_min_flux: %f units\n", lc_min_flux);
776 // printf("lc_flux_scaling_constant: %f 10^39 erg s^-1\n", lc_flux_scaling_constant);
777 // printf("lc_P_wind_time_scaling_constant: %f s\n", lc_P_wind_time_scaling_constant);
778 // printf("lc_ulx_lim: %f units \n", lc_ulx_lim);
779 // printf("lc_classification: %d \n", lc_classification);
780 //
781 // printf("\n");
782 //
783 // printf("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ \n");
784 // printf("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ \n");
785 // printf("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ \n");
786
787 // (%d, %f, %d, %d, %f, %f, %f, %d, ?)
788 // Insert CLASSIFICATIONS results into sqlite table
789 sql = sqlite3_mprintf("INSERT INTO CLASSIFICATIONS VALUES (%d, %f, %d, %d, %f, %f, %f, %d, %Q);", system_row_id[N], system_theta[N], system_dincl[N], system_inclination[N], lc_min_flux, lc_max_flux, lc_flux_scaling_constant, lc_classification, run_id);
790
791 /* Execute SQL statement */
792 rc = sqlite3_open(filename_database, &db);
793 rc = sqlite3_exec(db, sql, callback, 0, &zErrMsg);
794
795 if( rc != SQLITE_OK ){
796 fprintf(stderr, "SQL error: %s\n", zErrMsg);
797 sqlite3_free(zErrMsg);
798 } else {
799 // fprintf(stdout, "CLASSIFICATIONS Records created successfully\n");
800 }
801 sqlite3_close(db);
802 }
803
804 return 0;
805}
806