/*------------------------------------------------------------*/ /* */ /* Program to compute Dopplergrams using a MDI-like algorithm */ /* Author: S. Couvidat (based on code by J. Schou) */ /* Version 1.4 April 2008 */ /* */ /* uses a MDI-like algorithm with 6 tuning positions */ /* averages the velocities returned by 1st and 2nd Fourier */ /* coefficient, and by LCP and RCP */ /* the code also estimates the Fe I linewidth, linedepth, and */ /* the continuum intensity */ /* NB: this is a PRELIMINARY VERSION. Nothing has been decided*/ /* yet regarding the exact filtergram sequence, the way the */ /* data must be saved, the way the look-up tables are saved...*/ /* NB: OpenMP version not implemented yet, but will be done */ /* */ /*------------------------------------------------------------*/ #include #include #include #include #include //#include //OpenMP header (not used for now...) #include //GNU Scientific Library header #include #undef I //I is the complex number (0,1). We un-define to avoid confusion char *module_name = "velocity"; //name of the module #define kRecSetIn "input_series" //names of the arguments of the module #define kDSOut "dopplergram" #define kDSOut2 "magnetogram" #define kDSOut3 "linedepth" #define kDSOut4 "linewidth" #define kDSOut5 "continuum" //arguments of the module ModuleArgs_t module_args[] = { {ARG_STRING, kRecSetIn, "", "Input data series."}, {ARG_STRING, kDSOut, "", "Output Dopplergram series."}, {ARG_STRING, kDSOut2, "", "Output magnetogram series."}, {ARG_STRING, kDSOut3, "", "Output linedepth series."}, {ARG_STRING, kDSOut4, "", "Output linewidth series."}, {ARG_STRING, kDSOut5, "", "Output continuum intensity series."}, {ARG_END} }; /*-------------------------------------------------------------*/ /* DoIt is the entry point of the module */ /* the name MUST be DoIt for a DRMS module */ /*-------------------------------------------------------------*/ int DoIt(void) { char *inRecQuery = cmdparams_get_str(&cmdparams, kRecSetIn, NULL); //cmdparams is defined in jsoc_main.h char *dsout = cmdparams_get_str(&cmdparams, kDSOut, NULL); //series name of the output Dopplergram series char *dsout2 = cmdparams_get_str(&cmdparams, kDSOut2,NULL); //series name of the output magnetogram series char *dsout3 = cmdparams_get_str(&cmdparams, kDSOut3,NULL); //series name of the output linedepth series char *dsout4 = cmdparams_get_str(&cmdparams, kDSOut4,NULL); //series name of the output linewidth series char *dsout5 = cmdparams_get_str(&cmdparams, kDSOut5,NULL); //series name of the output continuum series int status = DRMS_SUCCESS; int error = 0; //int nthreads, tid, i, chunk; //for OpenMP printf("START!\n"); /*#pragma omp parallel private(nthreads, tid) { tid = omp_get_thread_num(); if (tid == 0) { nthreads = omp_get_num_threads(); printf("Number of threads = %d\n", nthreads); //Number of threads for OpenMP } }*/ /*drms_series_exists(drms_env, inRecQuery, &status); if (status == DRMS_ERROR_UNKNOWNSERIES) { printf("Input series %s doesn't exist.\n",inRecQuery); exit(EXIT_FAILURE); } if (status == DRMS_SUCCESS) { printf("Input series %s exists and is being retrieved.\n",inRecQuery); }*/ drms_series_exists(drms_env, dsout, &status); //check whether or not output series exists //drms_env is defined in jsoc_main.h if (status == DRMS_ERROR_UNKNOWNSERIES) { printf("Output Dopplergram series %s doesn't exist\n",dsout); //if the output series does not exit exit(EXIT_FAILURE); //we exit the program } if (status == DRMS_SUCCESS) { printf("Output Dopplergram series %s exists.\n",dsout); } drms_series_exists(drms_env, dsout2, &status); //check whether or not output series exists //drms_env is defined in jsoc_main.h if (status == DRMS_ERROR_UNKNOWNSERIES) { printf("Output magnetogram series %s doesn't exist\n",dsout2); //if the output series does not exit exit(EXIT_FAILURE); //we exit the program } if (status == DRMS_SUCCESS) { printf("Output magnetogram series %s exists.\n",dsout2); } drms_series_exists(drms_env, dsout3, &status); //check whether or not output series exists //drms_env is defined in jsoc_main.h if (status == DRMS_ERROR_UNKNOWNSERIES) { printf("Output linedepth series %s doesn't exist\n",dsout3); //if the output series does not exit exit(EXIT_FAILURE); //we exit the program } if (status == DRMS_SUCCESS) { printf("Output linedepth series %s exists.\n",dsout3); } drms_series_exists(drms_env, dsout4, &status); //check whether or not output series exists //drms_env is defined in jsoc_main.h if (status == DRMS_ERROR_UNKNOWNSERIES) { printf("Output linewidth series %s doesn't exist\n",dsout4); //if the output series does not exit exit(EXIT_FAILURE); //we exit the program } if (status == DRMS_SUCCESS) { printf("Output linewidth series %s exists.\n",dsout4); } drms_series_exists(drms_env, dsout5, &status); //check whether or not output series exists //drms_env is defined in jsoc_main.h if (status == DRMS_ERROR_UNKNOWNSERIES) { printf("Output continuum intensity series %s doesn't exist\n",dsout5); //if the output series does not exit exit(EXIT_FAILURE); //we exit the program } if (status == DRMS_SUCCESS) { printf("Output continuum intensity series %s exists.\n",dsout5); } //open the records from the input series DRMS_RecordSet_t *data = drms_open_records(drms_env,inRecQuery,&status); if (status == DRMS_SUCCESS && data != NULL) { int nRecs = data->n; //number of records in the input series printf("Number of filtergrams to be read= %d \n",nRecs); if (nRecs != 24) //this number is currently 24 (12 filtergrams with front camera, 12 with side camera. BUT WILL MOST LIKELY CHANGE) { printf("Problem: program requires 24 filtergrams to produce the observables \n"); drms_close_records(data,DRMS_FREE_RECORD); exit(EXIT_FAILURE); } //the following paragraph might change if the filtergram sequence is altered DRMS_Record_t *rec0 = data->records[0]; //first HMI filtergram in LCP DRMS_Record_t *rec1 = data->records[4]; //second filtergram in LCP DRMS_Record_t *rec2 = data->records[8]; //third filtergram in LCP DRMS_Record_t *rec3 = data->records[12]; //fourth filtergram in LCP DRMS_Record_t *rec4 = data->records[16]; //fifth filtergram in LCP DRMS_Record_t *rec5 = data->records[20]; //sixth filtergram in LCP DRMS_Record_t *rec6 = data->records[2]; //first HMI filtergram in RCP DRMS_Record_t *rec7 = data->records[6]; //second filtergram in RCP DRMS_Record_t *rec8 = data->records[10]; //third filtergram in RCP DRMS_Record_t *rec9 = data->records[14]; //fourth filtergram in RCP DRMS_Record_t *rec10 = data->records[18]; //fifth filtergram in RCP DRMS_Record_t *rec11 = data->records[22]; //sixth filtergram in RCP DRMS_Segment_t *segin = NULL; DRMS_Array_t *arrin0 = NULL; DRMS_Array_t *arrin1 = NULL; DRMS_Array_t *arrin2 = NULL; DRMS_Array_t *arrin3 = NULL; DRMS_Array_t *arrin4 = NULL; DRMS_Array_t *arrin5 = NULL; DRMS_Array_t *arrin6 = NULL; DRMS_Array_t *arrin7 = NULL; DRMS_Array_t *arrin8 = NULL; DRMS_Array_t *arrin9 = NULL; DRMS_Array_t *arrin10 = NULL; DRMS_Array_t *arrin11 = NULL; DRMS_RecordSet_t *dataout = NULL; DRMS_RecordSet_t *dataout2= NULL; DRMS_RecordSet_t *dataout3= NULL; DRMS_RecordSet_t *dataout4= NULL; DRMS_RecordSet_t *dataout5= NULL; DRMS_Record_t *recout = NULL; DRMS_Segment_t *segout = NULL; DRMS_Array_t *arrout = NULL; //array that will contain the Dopplergram DRMS_Array_t *arrout2 = NULL; //array that will contain the magnetogram DRMS_Array_t *arrout3 = NULL; //array that will contain the linedepth map DRMS_Array_t *arrout4 = NULL; //array that will contain the linewidth map DRMS_Array_t *arrout5 = NULL; //array that will contain the continuum intensity map const char *keyname = "FSN"; //1st prime key of input data (TEMPORARY: THE KEYWORDS TO USE HAVE NOT BEEN DECIDED YET) const char *keyname2 = "T_OBS"; //2nd prime key of input data const char *primekey = "FSN_START"; //1st prime key of output data const char *primekey2 = "T_START"; //2nd prime key of output data const char *key3 = "FSN_STOP"; //extra keyword for output data const char *key4 = "T_STOP"; //extra keyword for output data int keyvalue = 0; int keyvalue3 = 0; DRMS_Type_Value_t keyvalue2; DRMS_Type_Value_t keyvalue4; char *stime; TIME interntime; TIME interntime2; DRMS_Type_t type; DRMS_Type_t type2 = DRMS_TYPE_TIME; long nElem; const int N = 6; //number of filtergrams to produce a Dopplergram const double lam0 = 6173.3433; //solar Fe I line central wavelength in Angstrom double magnetic; //conversion factor from (v_RCP - v_LCP) to B double pv1 = 0.; double pv2 = 0.; double f1LCPc = 0.; double f1RCPc = 0.; double f1LCPs = 0.; double f1RCPs = 0.; double f2LCPc = 0.; double f2RCPc = 0.; double f2LCPs = 0.; double f2RCPs = 0.; double vLCP = 0.; double vRCP = 0.; double v2LCP = 0.; double v2RCP = 0.; int i = 0 ; //for loops int jjj = 0 ; int iii = 0 ; double temp0 = 0.; double temp = 0.; double temp2 = 0.; double temp3 = 0.; double tempbis = 0.; double temp2bis = 0.; double temp3bis = 0.; double meanL = 0.; double meanR = 0.; int axisout[2] = {4096,4096}; //size of the output arrays. Needs to be made variable??? magnetic = 1.0/(2.0*4.67e-5*0.000061733433*2.5*299792458.0); //Lande factor=2.5 for Fe I line //read the first filtergram segin = drms_segment_lookupnum(rec0, 0); arrin0 = drms_segment_read(segin, segin->info->type, &status); short *arrinL0 = arrin0->data ; type = arrin0->type; //data type (short most likely) nElem = arrin0->axis[0]*arrin0->axis[1]; //filtergram sizes (4096*4096 most likely) printf("dimensions %d %d\n",arrin0->axis[0],arrin0->axis[1]); if (arrin0->axis[0] != 4096 || arrin0->axis[1] != 4096 ) { printf("Input filtergrams are not 4096x4096 \n"); drms_close_records(data,DRMS_FREE_RECORD); exit(EXIT_FAILURE); } printf("data type: %s \n",drms_type2str(type)); //print the data type printf("total number of elements: %d \n",nElem); //print the total number of elements keyvalue = drms_getkey_int(rec0,keyname,&status); keyvalue3= drms_getkey_int(rec5,keyname,&status); keyvalue2= drms_getkey(rec0,keyname2,&type2,&status); printf("FSN value %d\n",keyvalue); //value of the FSN of the 1st filtergram keyvalue4= drms_getkey(rec5,keyname2,&type2,&status); interntime = keyvalue2.time_val; interntime2= keyvalue4.time_val; //sprint_time(stime,interntime,zone,-1); //not working or problem with my installation of the DRMS? printf("T_OBS value %f\n",interntime); type = DRMS_TYPE_SHORT; //type of the output data: SHORT INTEGER (Needs to be variable?) arrout = drms_array_create(type,2,axisout,NULL,&status); //create the array that will contain the Dopplergram arrout2 = drms_array_create(type,2,axisout,NULL,&status); //create the array that will contain the magnetogram arrout3 = drms_array_create(type,2,axisout,NULL,&status); //create the array that will contain the linedepth arrout4 = drms_array_create(type,2,axisout,NULL,&status); //create the array that will contain the linewidth arrout5 = drms_array_create(type,2,axisout,NULL,&status); //create the array that will contain the continuum intensity short *lam0g = arrout->data ; //Dopplergram short *B0g = arrout2->data; //magnetogram short *Idg = arrout3->data; //linedepth short *widthg = arrout4->data; //linewidth short *I0g = arrout5->data; //continuum memset(arrout->data, 0, drms_array_size(arrout) ); //fill the Dopplergram with 0 memset(arrout2->data, 0, drms_array_size(arrout2)); memset(arrout3->data, 0, drms_array_size(arrout3)); memset(arrout4->data, 0, drms_array_size(arrout4)); memset(arrout5->data, 0, drms_array_size(arrout5)); printf("number of elements of output: %d\n",drms_array_size(arrout)/sizeof(short)); //read all the filtergrams segin = NULL; segin = drms_segment_lookupnum(rec1, 0); arrin1 = drms_segment_read(segin, segin->info->type, &status); short *arrinL1 = arrin1->data ; segin = NULL; segin = drms_segment_lookupnum(rec2, 0); arrin2 = drms_segment_read(segin, segin->info->type, &status); short *arrinL2 = arrin2->data ; segin = NULL; segin = drms_segment_lookupnum(rec3, 0); arrin3 = drms_segment_read(segin, segin->info->type, &status); short *arrinL3 = arrin3->data ; segin = NULL; segin = drms_segment_lookupnum(rec4, 0); arrin4 = drms_segment_read(segin, segin->info->type, &status); short *arrinL4 = arrin4->data ; segin = NULL; segin = drms_segment_lookupnum(rec5, 0); arrin5 = drms_segment_read(segin, segin->info->type, &status); short *arrinL5 = arrin5->data ; segin = NULL; segin = drms_segment_lookupnum(rec6, 0); arrin6 = drms_segment_read(segin, segin->info->type, &status); short *arrinR0 = arrin6->data ; segin = NULL; segin = drms_segment_lookupnum(rec7, 0); arrin7 = drms_segment_read(segin, segin->info->type, &status); short *arrinR1 = arrin7->data ; segin = NULL; segin = drms_segment_lookupnum(rec8, 0); arrin8 = drms_segment_read(segin, segin->info->type, &status); short *arrinR2 = arrin8->data ; segin = NULL; segin = drms_segment_lookupnum(rec9, 0); arrin9 = drms_segment_read(segin, segin->info->type, &status); short *arrinR3 = arrin9->data ; segin = NULL; segin = drms_segment_lookupnum(rec10, 0); arrin10 = drms_segment_read(segin, segin->info->type, &status); short *arrinR4 = arrin10->data ; segin = NULL; segin = drms_segment_lookupnum(rec11, 0); arrin11 = drms_segment_read(segin, segin->info->type, &status); short *arrinR5 = arrin11->data ; double FSRNB = 0.172; //FSR Narrow-Band Michelson, in Angstroms (WILL CHANGE ONCE THE VALUE IS ACCURATELY MEASURED) double inttune = (N-1.)/2.; //number of tuning intervals (should be 5) double dtune = FSRNB/inttune; //wavelength separation between each tuning position double dv = 299792458.0/lam0; //conversion factor from wavelength to velocity double dvtune = dtune*dv; double tune[6] = {-2.5,-1.5,-0.5,0.5,1.5,2.5}; //tuning positions double angle[6]= {-2.5,-1.5,-0.5,0.5,1.5,2.5}; //double lutable[8]={-43.856043, 6639.640053, 3.525622, -149.660537, -8.698054, -35.306646, -1.301304, -13.029174}; //the look-up "table" for the first Fourier coefficient has been fitted by Legendre polynomials. error on v<0.85 m/s in dynamic range //double lutable2[16]={36.555202, 6893.883128, -25.139495, -43.734764, 15.977831, -70.980713, -12.871190, 50.768014, 1.559605, 51.709238, 13.825859, -119.781728, -14.057757, 77.764177, 6.177970, -25.871270}; //the look-up "table" for the second Fourier coefficient, fitted by Legendre polynomials. error v<0.85 m/s double cosi[6]; double sini[6]; double cos2i[6]; double sin2i[6]; double integralfilter0 = 0.049575435; //integral value of the HMI 1st filter (MIGHT CHANGE IF CALIBRATION PROVIDES NEW VALUES) double integralfilter1 = 0.049608823; //integral value of the HMI 1st filter (MIGHT CHANGE IF CALIBRATION PROVIDES NEW VALUES) double integralfilter2 = 0.049633797; //integral value of the HMI 1st filter (MIGHT CHANGE IF CALIBRATION PROVIDES NEW VALUES) double integralfilter3 = 0.049572778; //integral value of the HMI 1st filter (MIGHT CHANGE IF CALIBRATION PROVIDES NEW VALUES) double integralfilter4 = 0.049447747; //integral value of the HMI 1st filter (MIGHT CHANGE IF CALIBRATION PROVIDES NEW VALUES) double integralfilter5 = 0.049435261; //integral value of the HMI 1st filter (MIGHT CHANGE IF CALIBRATION PROVIDES NEW VALUES) double meannorm; meannorm=(integralfilter0+integralfilter1+integralfilter2+integralfilter3+integralfilter4+integralfilter5)/5.0; double period = (N-1)*dtune; double L0 = 0.0; double L1 = 0.0; double L2 = 0.0; double L3 = 0.0; double L4 = 0.0; double L5 = 0.0; double R0 = 0.0; double R1 = 0.0; double R2 = 0.0; double R3 = 0.0; double R4 = 0.0; double R5 = 0.0; pv1 = dvtune*(double)(N-1); //(dvtune*inttune*2.0) pv2 = pv1/2.; for(i=0;i<=5;++i) { tune[i] = tune[i]*dtune; angle[i]= angle[i]*2.0*M_PI/(double)N; cosi[i] = cos(angle[i]); sini[i] = sin(angle[i]); cos2i[i]= cos(2.0*angle[i]); sin2i[i]= sin(2.0*angle[i]); } //Read look-up tables (NEXT PARAGRAPH MIGHT CHANGE WHEN FORMAT OF THESE TABLES IS DECIDED) //Lookup table for 1st Fourier coefficient int ntest = 501; //depends on how the look-up table was produced double dvtest; double vtestt[ntest]; //double dlam = 1.0/1.75e3; //DO NOT CHANGE: must be the same as the code used to produce look-up tables //double dlamdv = 2.059205672212074294E-5;//DO NOT CHANGE double poly[ntest] ;//look-up "table" for 1st Fourier coefficient double poly2[ntest];//look-up "table" for 2nd Fourier coefficient int index_lo=0; int index_hi=ntest-1; int index = 0; int j=0; //double vtestmax; char *lookup = "lookup.txt"; //the look-up table is currently saved as a text file. Might change. const char *format="%13lf %13lf %13lf\n"; FILE *fp; fp = fopen(lookup,"r"); if(fp ==NULL) { printf("Error: can't open look-up tables file.\n"); return 1; } for(i=0;i<=ntest-1;++i) { fscanf(fp,format,&vtestt[i],&poly[i],&poly2[i]); //format of look-up table: vtestt is input velocities, poly and poly2 are output velocities //printf("%f %f %f\n",vtestt[i],poly[i],poly2[i]); } fclose(fp); //dvtest = dlam/dlamdv;//DO NOT CHANGE //for(i=0;i<=ntest-1;++i) vtestt[i] = dvtest*((double)i-((double)ntest-1.0)/2.0); //vtestmax=vtestt[ntest-1]; /*for(j=0;j<=ntest-1;++j) { poly[j] = lutable[0]+lutable[1]*vtestt[j]/vtestmax; for(i=2;i<=7;++i) { temp = vtestt[j]/vtestmax; poly[j] = poly[j]+lutable[i]*gsl_sf_legendre_Pl(i,temp);//Legendre polynomials } }*/ //Lookup table for 2nd Fourier coefficient /*for(j=0;j<=ntest-1;++j) { poly2[j] = lutable2[0]+lutable2[1]*vtestt[j]/vtestmax; for(i=2;i<=15;++i) { temp = vtestt[j]/vtestmax; poly2[j] = poly2[j]+lutable2[i]*gsl_sf_legendre_Pl(i,temp); } }*/ //Loop over all the pixels of the filtergrams do { do { if (arrinL0[iii*arrin0->axis[0]+jjj] != -32768) //-32768 is the value assigned to missing data. MIGHT CHANGE? { L0 = (double)arrinL0[iii*arrin0->axis[0]+jjj]/integralfilter0*meannorm; L1 = (double)arrinL1[iii*arrin0->axis[0]+jjj]/integralfilter1*meannorm; L2 = (double)arrinL2[iii*arrin0->axis[0]+jjj]/integralfilter2*meannorm; L3 = (double)arrinL3[iii*arrin0->axis[0]+jjj]/integralfilter3*meannorm; L4 = (double)arrinL4[iii*arrin0->axis[0]+jjj]/integralfilter4*meannorm; L5 = (double)arrinL5[iii*arrin0->axis[0]+jjj]/integralfilter5*meannorm; R0 = (double)arrinR0[iii*arrin0->axis[0]+jjj]/integralfilter0*meannorm; R1 = (double)arrinR1[iii*arrin0->axis[0]+jjj]/integralfilter1*meannorm; R2 = (double)arrinR2[iii*arrin0->axis[0]+jjj]/integralfilter2*meannorm; R3 = (double)arrinR3[iii*arrin0->axis[0]+jjj]/integralfilter3*meannorm; R4 = (double)arrinR4[iii*arrin0->axis[0]+jjj]/integralfilter4*meannorm; R5 = (double)arrinR5[iii*arrin0->axis[0]+jjj]/integralfilter5*meannorm; //First Fourier coefficients f1LCPc = cosi[0]*L0+cosi[1]*L1+cosi[2]*L2+cosi[3]*L3+cosi[4]*L4+cosi[5]*L5; f1RCPc = cosi[0]*R0+cosi[1]*R1+cosi[2]*R2+cosi[3]*R3+cosi[4]*R4+cosi[5]*R5; f1LCPs = sini[0]*L0+sini[1]*L1+sini[2]*L2+sini[3]*L3+sini[4]*L4+sini[5]*L5; f1RCPs = sini[0]*R0+sini[1]*R1+sini[2]*R2+sini[3]*R3+sini[4]*R4+sini[5]*R5; vLCP = atan2(-f1LCPs,-f1LCPc)*pv1/2.0/M_PI; //-f1LCPs and -f1LCPc because the first Fourier coefficient is -2/T*... vRCP = atan2(-f1RCPs,-f1RCPc)*pv1/2.0/M_PI; //Second Fourier coefficients f2LCPc = cos2i[0]*L0+cos2i[1]*L1+cos2i[2]*L2+cos2i[3]*L3+cos2i[4]*L4+cos2i[5]*L5; f2RCPc = cos2i[0]*R0+cos2i[1]*R1+cos2i[2]*R2+cos2i[3]*R3+cos2i[4]*R4+cos2i[5]*R5; f2LCPs = sin2i[0]*L0+sin2i[1]*L1+sin2i[2]*L2+sin2i[3]*L3+sin2i[4]*L4+sin2i[5]*L5; f2RCPs = sin2i[0]*R0+sin2i[1]*R1+sin2i[2]*R2+sin2i[3]*R3+sin2i[4]*R4+sin2i[5]*R5; v2LCP = atan2(-f2LCPs,-f2LCPc)*pv2/2.0/M_PI; v2RCP = atan2(-f2RCPs,-f2RCPc)*pv2/2.0/M_PI; v2LCP = fmod((v2LCP-vLCP+10.5*pv2),pv2)-pv2/2.0+vLCP; //we use the uncorrected velocity, i.e. phase, of the 1st Fourier coefficient to correct for the estimate of v2LCP and v2RCP, because the range of velocities obtained with the second Fourier coefficient is half the range of the first Fourier coefficient v2RCP = fmod((v2RCP-vRCP+10.5*pv2),pv2)-pv2/2.0+vRCP; //We linearly (might change) interpolate in the look-up table for the 1st Fourier coefficient to retrieve the actual velocities index = gsl_interp_bsearch(poly,vLCP,index_lo,index_hi); //returns the index of the array poly for which poly[index] <= vLCP <= poly[index+1] vLCP = vtestt[index]+(vLCP-poly[index])*(vtestt[index+1]-vtestt[index])/(poly[index+1]-poly[index]); index = gsl_interp_bsearch(poly,vRCP,index_lo,index_hi); vRCP = vtestt[index]+(vRCP-poly[index])*(vtestt[index+1]-vtestt[index])/(poly[index+1]-poly[index]); //We linearly (might change) interpolate in the look-up table for the 2nd Fourier coefficient to retrieve the actual velocities index = gsl_interp_bsearch(poly2,v2LCP,index_lo,index_hi); v2LCP = vtestt[index]+(v2LCP-poly2[index])*(vtestt[index+1]-vtestt[index])/(poly2[index+1]-poly2[index]); index = gsl_interp_bsearch(poly2,v2RCP,index_lo,index_hi); v2RCP = vtestt[index]+(v2RCP-poly2[index])*(vtestt[index+1]-vtestt[index])/(poly2[index+1]-poly2[index]); f1LCPc = -f1LCPc*dtune/period*2.0; f1RCPc = -f1RCPc*dtune/period*2.0; f1LCPs = -f1LCPs*dtune/period*2.0; f1RCPs = -f1RCPs*dtune/period*2.0; f2LCPc = -f2LCPc*dtune/period*2.0; f2RCPc = -f2RCPc*dtune/period*2.0; f2LCPs = -f2LCPs*dtune/period*2.0; f2RCPs = -f2RCPs*dtune/period*2.0; //We compute the Doppler velocity temp0 = (vLCP+vRCP+v2LCP+v2RCP)/4.; lam0g[iii*axisout[0]+jjj] = round(temp0); //simple average. Need weights? //We compute the l.o.s. magnetic field B0g[iii*axisout[0]+jjj] = round((vRCP-vLCP+v2RCP-v2LCP)/2.0*magnetic); //We compute the linewidth (in Angstroms) and multiply it by 10000 (NEED TO OPTIMIZE) because the type is SHORT INT temp = period/M_PI*sqrt(1.0/6.0*log((f1LCPc*f1LCPc+f1LCPs*f1LCPs)/(f2LCPc*f2LCPc+f2LCPs*f2LCPs))); tempbis = period/M_PI*sqrt(1.0/6.0*log((f1RCPc*f1RCPc+f1RCPs*f1RCPs)/(f2RCPc*f2RCPc+f2RCPs*f2RCPs))); widthg[iii*axisout[0]+jjj] = round(10000.*(temp+tempbis)*sqrt(log(2.0))); //we want the FWHM not the sigma of the Gaussian //We compute the linedepth temp2 = period/2.0*sqrt(f1LCPc*f1LCPc+f1LCPs*f1LCPs)/sqrt(M_PI)/temp*exp(M_PI*M_PI*temp*temp/period/period); temp2bis = period/2.0*sqrt(f1RCPc*f1RCPc+f1RCPs*f1RCPs)/sqrt(M_PI)/tempbis*exp(M_PI*M_PI*tempbis*tempbis/period/period); Idg[iii*axisout[0]+jjj] = round((temp2+temp2bis)/2.0); //We compute the continuum intensity temp3 = (vLCP+v2LCP)/2.0/dv; temp3bis = (vRCP+v2RCP)/2.0/dv; meanL = (L0+L1+L2+L3+L4+L5)/6.0; meanR = (R0+R1+R2+R3+R4+R5)/6.0; meanL = meanL + temp2/6.0*(exp(-(tune[0]-temp3)*(tune[0]-temp3)/temp/temp)+exp(-(tune[1]-temp3)*(tune[1]-temp3)/temp/temp)+exp(-(tune[2]-temp3)*(tune[2]-temp3)/temp/temp)+exp(-(tune[3]-temp3)*(tune[3]-temp3)/temp/temp)+exp(-(tune[4]-temp3)*(tune[4]-temp3)/temp/temp)+exp(-(tune[5]-temp3)*(tune[5]-temp3)/temp/temp)); meanR = meanR + temp2bis/6.0*(exp(-(tune[0]-temp3bis)*(tune[0]-temp3bis)/tempbis/tempbis)+exp(-(tune[1]-temp3bis)*(tune[1]-temp3bis)/tempbis/tempbis)+exp(-(tune[2]-temp3bis)*(tune[2]-temp3bis)/tempbis/tempbis)+exp(-(tune[3]-temp3bis)*(tune[3]-temp3bis)/tempbis/tempbis)+exp(-(tune[4]-temp3bis)*(tune[4]-temp3bis)/tempbis/tempbis)+exp(-(tune[5]-temp3bis)*(tune[5]-temp3bis)/tempbis/tempbis)); I0g[iii*axisout[0]+jjj] = round((meanL+meanR)/2.0); } } while (++jjj < axisout[0]); jjj = 0; } while (++iii < axisout[1]); drms_free_array(arrin0); drms_free_array(arrin1); drms_free_array(arrin2); drms_free_array(arrin3); drms_free_array(arrin4); drms_free_array(arrin5); drms_free_array(arrin6); drms_free_array(arrin7); drms_free_array(arrin8); drms_free_array(arrin9); drms_free_array(arrin10); drms_free_array(arrin11); //create a record in the series dsout (Dopplergrams) dataout = drms_create_records(drms_env,1,dsout ,DRMS_PERMANENT,&status); if (status != DRMS_SUCCESS) { printf("Could not create a record for the Dopplergram\n"); exit(EXIT_FAILURE); } if (status == DRMS_SUCCESS) { printf("Writing a record on the DRMS for the Dopplergram\n"); recout = dataout->records[0]; status = drms_setkey(recout,primekey,DRMS_TYPE_INT,&keyvalue); //create the prime key to identify the record. With my series su_couvidat.TestData, I decided that FSN_START is a prime key, along with T_START status = drms_setkey(recout,primekey2,DRMS_TYPE_TIME,&interntime); status = drms_setkey(recout,key3,DRMS_TYPE_INT,&keyvalue3); status = drms_setkey(recout,key4,DRMS_TYPE_TIME,&interntime2); printf("done\n"); segout = drms_segment_lookupnum(recout, 0); drms_segment_write(segout, arrout, 0); //write the file containing lam0g drms_close_records(dataout, DRMS_INSERT_RECORD); //insert the record in DRMS } //create a record in the series dsout2 (magnetograms) dataout2= drms_create_records(drms_env,1,dsout2,DRMS_PERMANENT,&status); if (status != DRMS_SUCCESS) { printf("Could not create a record for the magnetogram\n"); exit(EXIT_FAILURE); } if (status == DRMS_SUCCESS) { printf("Writing a record on the DRMS for the magnetogram\n"); recout = dataout2->records[0]; status = drms_setkey(recout,primekey,DRMS_TYPE_INT,&keyvalue); //create the prime key to identify the record. With my series su_couvidat.TestData2, I decided that FSN_START is a prime key, along with T_START status = drms_setkey(recout,primekey2,DRMS_TYPE_TIME,&interntime); status = drms_setkey(recout,key3,DRMS_TYPE_INT,&keyvalue3); status = drms_setkey(recout,key4,DRMS_TYPE_TIME,&interntime2); printf("done\n"); segout = drms_segment_lookupnum(recout, 0); drms_segment_write(segout, arrout2, 0); //write the file containing B0g drms_close_records(dataout2, DRMS_INSERT_RECORD); //insert the record in DRMS } //create a record in the series dsout3 (linedepths) dataout3= drms_create_records(drms_env,1,dsout3,DRMS_PERMANENT,&status); if (status != DRMS_SUCCESS) { printf("Could not create a record for the linedepth\n"); exit(EXIT_FAILURE); } if (status == DRMS_SUCCESS) { printf("Writing a record on the DRMS for the linedepth\n"); recout = dataout3->records[0]; status = drms_setkey(recout,primekey,DRMS_TYPE_INT,&keyvalue); //create the prime key to identify the record. With my series su_couvidat.TestData3, I decided that FSN_START is a prime key, along with T_START status = drms_setkey(recout,primekey2,DRMS_TYPE_TIME,&interntime); status = drms_setkey(recout,key3,DRMS_TYPE_INT,&keyvalue3); status = drms_setkey(recout,key4,DRMS_TYPE_TIME,&interntime2); printf("done\n"); segout = drms_segment_lookupnum(recout, 0); drms_segment_write(segout, arrout3, 0); //write the file containing the linedepth drms_close_records(dataout3, DRMS_INSERT_RECORD); //insert the record in DRMS } //create a record in the series dsout4 (linewidths) dataout4= drms_create_records(drms_env,1,dsout4,DRMS_PERMANENT,&status); if (status != DRMS_SUCCESS) { printf("Could not create a record for the linewidth\n"); exit(EXIT_FAILURE); } if (status == DRMS_SUCCESS) { printf("Writing a record on the DRMS for the linewidth\n"); recout = dataout4->records[0]; status = drms_setkey(recout,primekey,DRMS_TYPE_INT,&keyvalue); //create the prime key to identify the record. With my series su_couvidat.TestData4, I decided that FSN_START is a prime key, along with T_START status = drms_setkey(recout,primekey2,DRMS_TYPE_TIME,&interntime); status = drms_setkey(recout,key3,DRMS_TYPE_INT,&keyvalue3); status = drms_setkey(recout,key4,DRMS_TYPE_TIME,&interntime2); printf("done\n"); segout = drms_segment_lookupnum(recout, 0); drms_segment_write(segout, arrout4, 0); //write the file containing the linewidth drms_close_records(dataout4, DRMS_INSERT_RECORD); //insert the record in DRMS } //create a record in the series dsout5 (continuum intensity) dataout5= drms_create_records(drms_env,1,dsout5,DRMS_PERMANENT,&status); if (status != DRMS_SUCCESS) { printf("Could not create a record for the continuum intensity\n"); exit(EXIT_FAILURE); } if (status == DRMS_SUCCESS) { printf("Writing a record on the DRMS for the continuum intensity\n"); recout = dataout5->records[0]; status = drms_setkey(recout,primekey,DRMS_TYPE_INT,&keyvalue); //create the prime key to identify the record. With my series su_couvidat.TestData5, I decided that FSN_START is a prime key, along with T_START status = drms_setkey(recout,primekey2,DRMS_TYPE_TIME,&interntime); status = drms_setkey(recout,key3,DRMS_TYPE_INT,&keyvalue3); status = drms_setkey(recout,key4,DRMS_TYPE_TIME,&interntime2); printf("done\n"); segout = drms_segment_lookupnum(recout, 0); drms_segment_write(segout, arrout5, 0); //write the file containing the continuum intensity drms_close_records(dataout5, DRMS_INSERT_RECORD); //insert the record in DRMS } printf("COMPLETED!\n"); } if (data != NULL) { printf("close DRMS session \n"); drms_close_records(data,DRMS_FREE_RECORD); } return error; }