#include <stdio.h>
#include <math.h>
#include <stdlib.h>

main(argc,argv)
  int argc;
  char *argv[];
{
	double *time_bin, *count, *count_quad0, *count_quad1, *count_quad2, *count_quad3;
	double *count_4quads;
	double bg_num1_4quads,bg_num2_4quads,fg_num_4quads,mu_fg_4quads,epsilon_squ_4quads,signal_to_noise_4quads, signal_to_noise_4quads_lim;
	double time_num,time_num_all;
	double count_num;
	double signal_to_noise,signal_to_noise_lim,signal_to_noise_max,signal_to_noise_4quads_snrmax;
	double bg_avg,bg_num1,bg_num2;
	double fg_avg,fg_num;
	double epsilon_squ, mu_fg;
	double net_fg, image_0, image_1, image_2, image_3;
	double qd0, qd1, qd2, qd3;
	double trigger_time_step, trig_sig;
	double fg_dur_cal,bg_dur1_cal,bg_dur2_cal;
	double trigger_time_start, trigger_time_stop;
	double fg_count, T_90_start, T_90_stop, T90, frac;
	double bg_dur1_sec,fg_dur_sec,bg_dur2_sec,elapse_dur_sec, t_skip_sec, bin_size;
	int k1, T_90_start_bin, T_90_stop_bin;
	int trigger_time_start_bin, trigger_time_stop_bin;
	int t_bg2_start,t_bg2_stop;
	int t_fg_start,t_fg_stop;
	int bg_dur1,bg_dur2,fg_dur,elapse_dur;
	int i,imax,j,line_count,size,t_step,t_skip;
	int time_range, t_final,i_real;
	int energy_i,energy_num;
	int line_i, real_size;	
	int k,zero_check;
	int i_quad0, i_quad1, i_quad2, i_quad3;
	int trigger_mode_flag, rate_threshold_flag, image_threshold_flag;
	char inputname_quad0[500],inputname_quad1[500],inputname_quad2[500],inputname_quad3[500];
	char home[1000],foldername[500];
	char inputname_short[500];
	char line[5000],criterio[1000];
	FILE *input_quad0,*input_quad1,*input_quad2,*input_quad3;

	if (argc!=18) {
          printf("simple_trigger <inputfile> <bg_dur1(bin)> <fg_dur(bin)> <bg_dur2(bin)> <elapse_dur(bin)> <t_skip(bin)> <signal_to_noise_lim_squ> <qd0> <qd1> <qd2> <qd3> <energy_num> <criterio> <lc_folder> <bin_size> <image_threshold> <trigger_mode_flag>\n");
          return 0;
        }

	sscanf(argv[1],"%s",inputname_short);
	sscanf(argv[2],"%lf",&bg_dur1_sec);
	sscanf(argv[3],"%lf",&fg_dur_sec);
	sscanf(argv[4],"%lf",&bg_dur2_sec);
	sscanf(argv[5],"%lf",&elapse_dur_sec);
	sscanf(argv[6],"%lf",&t_skip_sec);
	sscanf(argv[7],"%lf",&trig_sig);
	signal_to_noise_lim = sqrt(trig_sig);
	sscanf(argv[8],"%lf",&qd0);
	sscanf(argv[9],"%lf",&qd1);
	sscanf(argv[10],"%lf",&qd2);
	sscanf(argv[11],"%lf",&qd3);
	sscanf(argv[12],"%d",&energy_num);
	sscanf(argv[13],"%s",criterio);
	sscanf(argv[14],"%s",home);
	sscanf(argv[15],"%lf",&bin_size);
	sscanf(argv[16],"%lf",&signal_to_noise_4quads_lim);
	sscanf(argv[17],"%d",&trigger_mode_flag);

	//change duration unit from sec to bin_size
        bg_dur1 = round(bg_dur1_sec/bin_size);
        fg_dur = round(fg_dur_sec/bin_size);
        bg_dur2 = round(bg_dur2_sec/bin_size);
        elapse_dur = round(elapse_dur_sec/bin_size);
        t_skip = round(t_skip_sec/bin_size);

	//print out inputs
	printf("#########################\n");
	printf("## Running trigger code with the following settings:\n");
	printf("## Light curve file: %s\n",inputname_short);
	if(trigger_mode_flag==1){
		printf("## Trigger mode: Long trigger\n");
	}
	else if(trigger_mode_flag==0){
                printf("## Trigger mode: Short trigger\n");
        }
	else{printf("Wrong! Trigger mode does not exist. End process.\n"); return 0;}
	if(energy_num==0){
		printf("## Energy band: 15-25 keV\n");
	}
	if(energy_num==1){
                printf("## Energy band: 15-50 keV\n");
        }
	if(energy_num==2){
                printf("## Energy band: 25-100 keV\n");
        }
	if(energy_num==3){
                printf("## Energy band: 50_350 keV\n");
        }	
	printf("## Time bracket: bg1_dur=%d, fg_dur=%d, bg2_dur=%d, elapse_dur=%d (unit: sec)\n",bg_dur1,fg_dur,bg_dur2,elapse_dur);
	printf("## Quadrant info: qd0=%d, qd1=%d, qd2=%d, qd3=%d\n",(int)qd0,(int)qd1,(int)qd2,(int)qd3);
	printf("## Signal-to-noise ratio threshold = %lf\n",signal_to_noise_lim);
	printf("## Image threshold = %lf\n",signal_to_noise_4quads_lim);
	printf("#########################\n");

	//read in light curve filenames of the input energy band
	if(energy_num == 0){
        sprintf(inputname_quad0,"%s/%s_bat_15_25kev_lc_sim_quad0.txt",home,inputname_short);
	sprintf(inputname_quad1,"%s/%s_bat_15_25kev_lc_sim_quad1.txt",home,inputname_short);
        sprintf(inputname_quad2,"%s/%s_bat_15_25kev_lc_sim_quad2.txt",home,inputname_short);
	sprintf(inputname_quad3,"%s/%s_bat_15_25kev_lc_sim_quad3.txt",home,inputname_short);
	}
	else if(energy_num == 1){
        sprintf(inputname_quad0,"%s/%s_bat_15_50kev_lc_sim_quad0.txt",home,inputname_short);
	sprintf(inputname_quad1,"%s/%s_bat_15_50kev_lc_sim_quad1.txt",home,inputname_short);
        sprintf(inputname_quad2,"%s/%s_bat_15_50kev_lc_sim_quad2.txt",home,inputname_short);
	sprintf(inputname_quad3,"%s/%s_bat_15_50kev_lc_sim_quad3.txt",home,inputname_short);
	}
	else if(energy_num == 2){
        sprintf(inputname_quad0,"%s/%s_bat_25_100kev_lc_sim_quad0.txt",home,inputname_short);
	sprintf(inputname_quad1,"%s/%s_bat_25_100kev_lc_sim_quad1.txt",home,inputname_short);
	sprintf(inputname_quad2,"%s/%s_bat_25_100kev_lc_sim_quad2.txt",home,inputname_short);
	sprintf(inputname_quad3,"%s/%s_bat_25_100kev_lc_sim_quad3.txt",home,inputname_short);
	}
	else if(energy_num == 3){
	sprintf(inputname_quad0,"%s/%s_bat_50_350kev_lc_sim_quad0.txt",home,inputname_short);
	sprintf(inputname_quad1,"%s/%s_bat_50_350kev_lc_sim_quad1.txt",home,inputname_short);
	sprintf(inputname_quad2,"%s/%s_bat_50_350kev_lc_sim_quad2.txt",home,inputname_short);
	sprintf(inputname_quad3,"%s/%s_bat_50_350kev_lc_sim_quad3.txt",home,inputname_short);
	}
	else{
	printf("Wrong! Energy_num = %d, need to be either 0 ,1, 2, or 3. End process.\n",energy_num); return 0;
	}

	//count the number of line in inputfile
	input_quad0=fopen(inputname_quad0,"r");
	line_count=0;
	while(fgets(line,sizeof(line),input_quad0)!=NULL){
               if(strncmp(line,"#",1)!=0){
			line_count++;
 	       }//if strncmp
        }//while
	fclose(input_quad0);


	//read in file
	input_quad0=fopen(inputname_quad0,"r");
        input_quad1=fopen(inputname_quad1,"r");
        input_quad2=fopen(inputname_quad2,"r");
        input_quad3=fopen(inputname_quad3,"r");
	size=line_count;
	time_bin=(double *)calloc(size,sizeof(double));
 	count=(double *)calloc(size,sizeof(double));
	count_quad0=(double *)calloc(size,sizeof(double));
	count_quad1=(double *)calloc(size,sizeof(double));
	count_quad2=(double *)calloc(size,sizeof(double));
	count_quad3=(double *)calloc(size,sizeof(double));
	count_4quads=(double *)calloc(size,sizeof(double));

	//get data from light curve files
	i_quad0=0;
	while(fgets(line,sizeof(line),input_quad0)!=NULL){       
		if(strncmp(line,"#",1)!=0){ 
			sscanf(line,"%lf %lf\n",&time_num,&count_num);
			time_bin[i_quad0]=time_num;
			count_quad0[i_quad0] = count_num;
			i_quad0++;
		}//if strncmp
	}//while
	i_quad1=0;
	while(fgets(line,sizeof(line),input_quad1)!=NULL){
                if(strncmp(line,"#",1)!=0){
                        sscanf(line,"%lf %lf\n",&time_num,&count_num);
                        count_quad1[i_quad1] = count_num;
                        i_quad1++;
                }//if strncmp
        }//while	
	i_quad2=0;
        while(fgets(line,sizeof(line),input_quad2)!=NULL){
                if(strncmp(line,"#",1)!=0){
                        sscanf(line,"%lf %lf\n",&time_num,&count_num);
                        count_quad2[i_quad2] = count_num;
                        i_quad2++;
                }//if strncmp
        }//while  
	i_quad3=0;
        while(fgets(line,sizeof(line),input_quad3)!=NULL){
                if(strncmp(line,"#",1)!=0){
                        sscanf(line,"%lf %lf\n",&time_num,&count_num);
                        count_quad3[i_quad3] = count_num;
                        i_quad3++;
                }//if strncmp
        }//while  

	//add up quadrant light curves based on input quadrant info
	for(i=0;i<size;i++){
		count[i]=qd0*count_quad0[i]+qd1*count_quad1[i]+qd2*count_quad2[i]+qd3*count_quad3[i];
	}

	//get the 4quads count for image trigger
	i_real = 0;
	for(i=0;i<size;i++){
		i_real += 1;
        	count_4quads[i]=count_quad0[i]+count_quad1[i]+count_quad2[i]+count_quad3[i];
        }

	//give warning messages if elapse_dur = 0
	if(elapse_dur == 0){printf("Warning: Elapse time is zero!\n");}


	//check if fg_dur ==0, if so, end calcultation
	if(fg_dur == 0){
                printf("Wrong! duration for fg cannot be zero! End process.\n"); return 0;
        }

	printf("## Results:\n");
	printf("## criterion signal_to_noise_ratio fg_start_time fg_end_time bg_start_time bg_end_time bg2_start_time bg2_end_time signal_to_noise_ratio_4quads\n");	

	//if fg_dur!=0, start calculation
	if(fg_dur!=0){
		signal_to_noise_max = 0.0;
		signal_to_noise_4quads_snrmax = 0.0;
		rate_threshold_flag = 0;
		image_threshold_flag = 0;
		real_size = i_real;
		time_range = bg_dur1 + elapse_dur + fg_dur + elapse_dur + bg_dur2;	
	
		//make sure the light curve is long enough for the time_range
		if((real_size-time_range+1)<0){printf("Wrong! Light curve is shorter than the trigger time bracket! End process.\n"); return 0;}

		//go through the light curve to find triggers
		for(t_step=0;t_step<(real_size-time_range+1);t_step+=t_skip){
			t_final = t_step+time_range;
		
			//see if there are any zero or missing times in the range, if so, skip the time bracket
			zero_check = 1;
			if((t_step-2)>=0 && (t_final+2)<real_size){
				for(k=(t_step-1);k<=(t_final+2);k++){
					if(count[k]==0.0){zero_check = 0; printf("### Wrong! Zero at bin %d = time %lf. Skip this time step.\n ", k, time_bin[k]); break;}
					else if((time_bin[k]-time_bin[k-1])>2.0*(time_bin[2]-time_bin[1])){zero_check = 0; printf("### Wrong! Time missing at bin %d, time %lf %lf. Skip this time step.\n", k, time_bin[k], time_bin[k-1]); break;} //check if there are missing times
				}//for
			}//if

			//calculate the signal-to-noise ratio in the time bracket 
			if((t_step+time_range)<=real_size && zero_check == 1){
				//count in bg1
				bg_num1=0.0;
				bg_num1_4quads=0.0;
				for(j=t_step;j<(t_step+bg_dur1);j++){
					bg_num1+=count[j];
					bg_num1_4quads+=count_4quads[j];
				}//for bg
				//count in bg2
				bg_num2=0.0;
				bg_num2_4quads=0.0;
				t_bg2_start = t_step+bg_dur1+elapse_dur+fg_dur+elapse_dur;
				t_bg2_stop = t_bg2_start+bg_dur2;
				for(j=t_bg2_start;j<t_bg2_stop;j++){
				        bg_num2+=count[j];
					bg_num2_4quads+=count_4quads[j];
                	       	 }//for bg
				
				//count in fg
				t_fg_start = t_step+bg_dur1+elapse_dur;
				t_fg_stop = t_fg_start+fg_dur;
				fg_num=0.0;
				fg_num_4quads=0.0;
				for(k=t_fg_start;k<t_fg_stop;k++){
					fg_num+=count[k];
					fg_num_4quads+=count_4quads[k];
				}//for fg


				//calculate signal-to-noise ratio based on different situations:
				//Long trigger:
				if(trigger_mode_flag == 1){
					//In the case of long trigger, 
					//SNR is calculated based on Carlo Graziani's ``Spiffy'' Trigger paper, eq 3.
					//long trigger; 2-sided background	
					if(bg_dur1!=0 && bg_dur2!=0){
						if(bg_num1!=0.0 && bg_num2!=0.0){
							//calculate snr based on Carlo's paper
							fg_dur_cal=(double)fg_dur;
							bg_dur1_cal=(double)bg_dur1;
							bg_dur2_cal=(double)bg_dur2;
							epsilon_squ = fg_dur_cal*fg_dur_cal/(bg_dur1_cal*bg_dur1_cal/bg_num1+bg_dur2_cal*bg_dur2_cal/bg_num2); 
							mu_fg = ((bg_dur1_cal+bg_dur2_cal)/fg_dur_cal)*epsilon_squ;
							signal_to_noise = (fg_num - mu_fg)/sqrt(mu_fg + epsilon_squ);
	
							//calculate 4quads number for image trigger
							epsilon_squ_4quads = fg_dur_cal*fg_dur_cal/(bg_dur1_cal*bg_dur1_cal/bg_num1_4quads+bg_dur2_cal*bg_dur2_cal/bg_num2_4quads);
							mu_fg_4quads = ((bg_dur1_cal+bg_dur2_cal)/fg_dur_cal)*epsilon_squ_4quads;				
							signal_to_noise_4quads = (fg_num_4quads - mu_fg_4quads)/sqrt(mu_fg_4quads + epsilon_squ_4quads);	
						}//if bg_num1!=0 and bg_num2!=0
						else{
							printf("Warning: either photon count in bg1 or bg2 is zero!!!; bg_num1=%e, bg_num2=%e\n",bg_num1, bg_num2);
						}//else if bg_num1==0 or bg_num2==0
					}//if bg_dur!=0
	
					//long trigger, 1-sided background, bg2=0
					if(bg_dur1!=0 && bg_dur2 == 0){
						if(bg_num1!=0.0){
							//this is the case of one-sided background (bg2=0)
      		                       			//calculate snr based on Carlo's paper, but drop all the terms related to bg_dur2
       	                        			fg_dur_cal=(double)fg_dur;
                                			bg_dur1_cal=(double)bg_dur1;
                                			epsilon_squ = fg_dur_cal*fg_dur_cal/(bg_dur1_cal*bg_dur1_cal/bg_num1);
                                			mu_fg = ((bg_dur1_cal)/fg_dur_cal)*epsilon_squ;
        	                        		signal_to_noise = (fg_num - mu_fg)/sqrt(mu_fg + epsilon_squ);
	
                		                	//calculate 4quads number for image trigger
                                			epsilon_squ_4quads = fg_dur_cal*fg_dur_cal/(bg_dur1_cal*bg_dur1_cal/bg_num1_4quads);
                                			mu_fg_4quads = ((bg_dur1_cal)/fg_dur_cal)*epsilon_squ_4quads;
                                			signal_to_noise_4quads = (fg_num_4quads - mu_fg_4quads)/sqrt(mu_fg_4quads + epsilon_squ_4quads);
						}//if bg_num1!=0
						else{
							printf("Warning: photon count in bg1 is zero!!!; bg_num1=%e\n",bg_num1);
						} 
                        		}//if bg_dur2=0

					//long trigger, 1-sided background, bg1=0
					if(bg_dur1==0 && bg_dur2 != 0){
                                        	if(bg_num2!=0.0){
							//this is the case of one-sided background (bg1=0)
                                        		//calculate snr based on Carlo's paper, but drop all the terms related to bg_dur1
							fg_dur_cal=(double)fg_dur;
                                        		bg_dur2_cal=(double)bg_dur2;
                                        		epsilon_squ = fg_dur_cal*fg_dur_cal/(bg_dur2_cal*bg_dur2_cal/bg_num2);
                                        		mu_fg = ((bg_dur2_cal)/fg_dur_cal)*epsilon_squ;
                                        		signal_to_noise = (fg_num - mu_fg)/sqrt(mu_fg + epsilon_squ);

							//calculate 4quads number for image trigger
							epsilon_squ_4quads = fg_dur_cal*fg_dur_cal/(bg_dur2_cal*bg_dur2_cal/bg_num2_4quads);
                                        		mu_fg_4quads = ((bg_dur2_cal)/fg_dur_cal)*epsilon_squ_4quads;
                                        		signal_to_noise_4quads = (fg_num_4quads - mu_fg_4quads)/sqrt(mu_fg_4quads + epsilon_squ_4quads);
                                        		}//if bg_num2!=0
						else{
                                                	printf("Warning: photon count in bg2 is zero!!!; bg_num2=%e\n",bg_num2);
                                        	}
                                	}//if bg_dur1=0					
				}//if long_trigger (trigger_mode_flag==1)				

				//short trigger:
				if(trigger_mode_flag==0){
					//this is the case of short trigger, which only has one background
					//calculate snr based on Fenimore et al. 2004 (eq. 1), but without the sigma_min term
					//This is similar to the long trigger and Carlo's equation with one background, but only has sigma_foreground term,
					//and does not has the sigma_background term
					if(bg_dur1!=0 && bg_dur2==0){
                               			if(bg_num1!=0.0){
							fg_dur_cal=(double)fg_dur;
	                       			        bg_dur1_cal=(double)bg_dur1;
                                			mu_fg = fg_dur_cal/bg_dur1_cal*bg_num1;
                                			signal_to_noise = (fg_num - mu_fg)/sqrt(mu_fg);

                                			//calculate 4quads number for image trigger
                                			mu_fg_4quads = fg_dur_cal/bg_dur1_cal*bg_num1_4quads;
                                			signal_to_noise_4quads = (fg_num_4quads - mu_fg_4quads)/sqrt(mu_fg_4quads);
                                			}//if bg_num1!=0
							else{
                                			        printf("Warning: photon count in bg1 is zero!!!; bg_num1=%e\n",bg_num1);
                                			}
                        		}//if bg_dur2=0
					if(bg_dur2!=0 && bg_dur1==0){
                                                if(bg_num2!=0.0){
                                                        fg_dur_cal=(double)fg_dur;
                                                        bg_dur2_cal=(double)bg_dur2;
                                                        mu_fg = fg_dur_cal/bg_dur2_cal*bg_num2;
                                                        signal_to_noise = (fg_num - mu_fg)/sqrt(mu_fg);
							
							//calculate 4quads number for image trigger
							mu_fg_4quads = fg_dur_cal/bg_dur2_cal*bg_num2_4quads;
                                                        signal_to_noise_4quads = (fg_num_4quads - mu_fg_4quads)/sqrt(mu_fg_4quads);
                                                        }//if bg_num2!=0
                                                        else{
                                                                printf("Warning: photon count in bg1 is zero!!!; bg_num2=%e\n",bg_num2);
                                                        }
                                        }//if bg_dur1=0
					if(bg_dur1!=0 && bg_dur2!=0){
						printf("Wrong! This is short trigger, either bg1 or bg2 has to be zero. End process.\n"); return 0;
					}

				}//if short_trigger (trigger_mode_flat==0)

				//get signal-to-noise above the input limit
				if(signal_to_noise > signal_to_noise_lim){
					rate_threshold_flag = 1;	
	
					//get signal-to-noise-4quads above the image threshold
					if(signal_to_noise_4quads > signal_to_noise_4quads_lim){
						image_threshold_flag = 1;
						printf("%s %e %e %e %e %e %e %e %e\n", criterio,signal_to_noise,time_bin[t_fg_start],time_bin[t_fg_stop],time_bin[t_step],time_bin[t_step+bg_dur1],time_bin[t_bg2_start],time_bin[t_bg2_stop], signal_to_noise_4quads);
	
						//get the max signal_to_noise
						if(signal_to_noise > signal_to_noise_max){
		                                        signal_to_noise_max = signal_to_noise;
       			                                trigger_time_step = time_bin[t_fg_start];
                                        		trigger_time_start = time_bin[t_step];
                                        		trigger_time_stop = time_bin[t_final];
                                        		trigger_time_start_bin = t_step;
                                        		trigger_time_stop_bin = t_final;
                                        		signal_to_noise_4quads_snrmax = signal_to_noise_4quads;
                                		}//if signal-to-noise-max

					} // if signal_to_noise_ratio_4quads

				}//if signal_to_noise_ratio
	
			}//if in range
	
		}//for t_step		
	}//if fg_dur!=0

	// see if the burst is above the required signal-to-noise and print out summary
	printf("## Summary of trigger criterion '%s':\n",criterio);
	if(rate_threshold_flag==1){
                if(image_threshold_flag==1){
                        printf("## triggered! Max signal-to-noise = %lf\n", signal_to_noise_max);
                        printf("## trigger time: %lf\n", trigger_time_step);
                        printf("## trigger bracket time range: from %lf to %lf\n", trigger_time_start, trigger_time_stop);
                        printf("#########################\n");
                }//if signal_to_noise_4quads
                else{
                        printf("## Not triggered: Pass rate trigger but do not pass the image threshold\n");
                        printf("#########################\n");
                }
        }//if
        else{printf("## Not triggered...\n");
                printf("#########################\n");
        }//else

	
	fclose(input_quad0);
	fclose(input_quad1);
	fclose(input_quad2);
	fclose(input_quad3);
	
}