JSOC_Documentation: proj/lev1.5_hmi/apps/HMI

00001  /*----------------------------------------------------------------------------------------------------------*/
00002 /* ANSI C99 CODE TO PRODUCE LOOK-UP TABLES FOR HMI                                                          */
00003 /* FOR THE OBSMODE ( BECAUSE WE TAKE INTO ACCOUNT THE SOLAR LINE PROFILE IN OBSMODE, WHICH IS DIFFERENT     */
00004 /* FROM CALMODE)                                                                                            */
00005 /* THEREFORE THE LOOK-UP TABLES ARE TECHNICALLY VALID ONLY FOR A GIVEN VALUE OF THE SOLAR RADIUS            */
00006 /* THIS VALUE IS THE solarradiustable VALUE OF HMIparam.h                                                   */
00007 /*                                                                                                          */
00008 /* Author: Sebastien Couvidat                                                                               */
00009 /* Mostly based on code by Jesper Schou                                                                     */
00010 /* version 3.0 MARCH 2011                                                                                   */
00011 /* version 3.1 JANUARY 2012                                                                                 */
00012 /* version 3.2 DECEMBER 2013: support for larger crop radii                                                 */
00013 /* version 3.3 January 2014: added calibration 13 effective at retune January 15, 2014                      */
00014 /*                                                                                                          */
00015 /* ASSUMPTION FOR THE LINE PROFILES: they are linear in theta, the angle from disk center                   */
00016 /*                                                                                                          */
00017 /* plus, run command line: limit stack unlimited                                                            */
00018 /* before running the code                                                                                  */
00019 /*                                                                                                          */
00020 /* cotune sequence from cotune_front_550752.txt (2007/09/05)                                                */
00021 /* for testing purpose: lookup phasemap="su_couvidat.phasemaps[707104]" lookup="su_couvidat.lookup" HCME1=37 HCMWB=59 HCMPOL=0 HCMNB=82 NUM=6 hcamid=1*/
00022 /* for testing purpose: lookup phasemap="su_couvidat.phasemaps[744586]" lookup="su_couvidat.lookup" HCME1=37 HCMWB=59 HCMPOL=0 HCMNB=82 NUM=6 hcamid=1*/
00023 /*                                                                                                          */
00024 /*                                                                                                          */
00025 /* CONVENTION FOR THE VELOCITIES:                                                                           */
00026 /* POSITIVE VELOCITIES CORRESPOND TO REDSHIFT (MOVEMENTS AWAY FROM THE OBSERVER)                            */
00027 /*                                                                                                          */
00028 /* NB: in this code the solar line is based on R. K. Ulrich line profile, while in the phasemaps.c code I   */
00029 /* use a Gaussian profile. That means an error on the central wavelength by about 2 mA, or about 97 m/s     */
00030 /* of systematic velocity error                                                                             */
00031 /* WE COULD PARTLY CORRECT FOR THIS ERROR BY SHIFTING ALL THE TEST VELOCITIES BY 97 m/s                     */
00032 /*                                                                                                          */
00033 /* NB: THIS CODE IS THE DRMS MODULE COUNTERPART OF THE LOOKUP2.PRO IDL CODE                                 */
00034 /* THE LOOK-UP TABLES RETURNED BY LOOKUP2.PRO DIFFER BY, AT MOST, 0.04 m/s                                  */
00035 /*                                                                                                          */
00036 /*----------------------------------------------------------------------------------------------------------*/
00037 
00038 /*
00039  * Changes for non-normal SDO orientation marked with comment including "ROLL 2018"
00040  * Need to read and write keywords CROTA2 and OBS_VR
00041  */
00042 
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include <math.h>
00046 #include <omp.h>
00047 #include <HMIparam.h>  //contains the FSRs of the HMI filter elements
00048 
00049 char *module_name    = "lookup";      //name of the module
00050 
00051 #define kRecSetIn      "phasemap" //names of the arguments of the module
00052 #define kDSOut         "lookup"
00053 #define HCM_E1         "HCME1"
00054 #define HCM_WB         "HCMWB"
00055 #define HCM_NB         "HCMNB"
00056 #define HCM_POL        "HCMPOL"
00057 #define NUM            "NUM"
00058 #define khcamid        "hcamid"       //front or side camera?                
00059 
00060 //convention for light and dark frames for keyword HCAMID
00061 #define LIGHT_SIDE  2   //SIDE CAMERA
00062 #define LIGHT_FRONT 3   //FRONT CAMERA
00063 #define DARK_SIDE   0   //SIDE CAMERA
00064 #define DARK_FRONT  1   //FRONT CAMERA
00065               
00066                         //arguments of the module
00067 ModuleArgs_t module_args[] =        
00068 {
00069      {ARG_STRING, kRecSetIn, "",  "Input phase-map series."},
00070      {ARG_STRING, kDSOut,    "",  "Lookup table series."},
00071      {ARG_INT   , NUM,       "",  "number of co-tuning positions"},
00072      {ARG_INT   , HCM_E1,    "",  "HCM position for E1 (in step units, should be at index 10 of cotune table)"},
00073      {ARG_INT   , HCM_WB,    "",  "HCM position for WB (in step units, should be at index 10 of cotune table)"},
00074      {ARG_INT   , HCM_NB,    "",  "HCM position for NB (in step units, should be at index 10 of cotune table)"},
00075      {ARG_INT   , HCM_POL,   "",  "HCM position for the tuning polarizer (in step units, should be at index 10 of cotune table)"},
00076      {ARG_INT   , khcamid,  "1",  "Front (1) or Side (0) camera?"},
00077      {ARG_INT   , "cal",      0,  "calibration used"},
00078      {ARG_END}
00079 };
00080 
00081 
00082 /*-----------------------------------------------------------------------------*/
00083 /* Function to perform linear interpolation                                    */
00084 /* found on the internet                                                       */
00085 /* returns the values yinterp at points x of 1D function yv (at the points xv) */
00086 /*-----------------------------------------------------------------------------*/
00087 
00088 void lininterp1f(double *yinterp, double *xv, double *yv, double *x, double ydefault, int m, int minterp)
00089 {
00090     int i, j; 
00091     int nrowsinterp, nrowsdata;
00092     nrowsinterp = minterp;
00093     nrowsdata = m;
00094     for (i=0; i<nrowsinterp; i++)
00095     {
00096             if((x[i] < xv[0]) || (x[i] > xv[nrowsdata-1]))
00097                 yinterp[i] = ydefault;
00098             else
00099             {   for(j=1; j<nrowsdata; j++)
00100             {      if(x[i]<=xv[j])
00101             {          yinterp[i] = (x[i]-xv[j-1]) / (xv[j]-xv[j-1]) * (yv[j]-yv[j-1]) + yv[j-1];
00102                        break;
00103             }
00104             }
00105             }
00106     }
00107 } 
00108 
00109 //to perform a basic integration using a piecewise linear scheme
00110 double integration(double *x,double *y,int N)
00111 {
00112   int i;
00113   double total=0.0;
00114 
00115   for(i=0;i<=N-2;++i) total=total+(x[i+1]-x[i])*(y[i+1]+y[i])/2.;
00116 
00117   return total;
00118               
00119 }
00120 
00121 //to perform a basic integration using the composite Simpson's rule
00122 //N must be odd so that the number of intervals (N-1) is even
00123 double simpson(double *x,double *y,int N)
00124 {
00125   int i;
00126   double h;
00127   double total=0.0;
00128 
00129   h=x[1]-x[0];
00130   total=h/3.0*(y[0]+y[N-1]);
00131   for(i=1;i<=(N-1)/2  ;++i) total+=h/3.0*4.0*y[2*i-1];
00132   for(i=1;i<=(N-1)/2-1;++i) total+=h/3.0*2.0*y[2*i];
00133 
00134   return total;
00135 } 
00136 
00137 
00138 /*---------------------------------------------------------------------------------------------------------------------------*/
00139 /*                                                                                                                           */
00140 /*  MAIN PROGRAM                                                                                                             */
00141 /*                                                                                                                           */
00142 /*---------------------------------------------------------------------------------------------------------------------------*/
00143 
00144 
00145 int DoIt(void) {
00146 
00147   int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
00148   int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
00149 
00150   int nthreads, tid, procs, maxt, inpar, dynamic, nested;//for OpenMP parallelization
00151   
00152   nthreads=omp_get_num_procs();                          //number of threads supported by the machine where the code is running
00153   omp_set_num_threads(nthreads);                         //set the number of threads to the maximum value
00154   printf("number of threads= %d\n",nthreads);
00155   
00156   int HCME1       = cmdparams_get_int(&cmdparams,HCM_E1, NULL); //tuning position of E1 correponding to center of Fe I line at rest (target wavelength)
00157   int HCMWB       = cmdparams_get_int(&cmdparams,HCM_WB, NULL); //tuning position of WB Michelson correponding to center of Fe I line at rest (target wavelength)
00158   int HCMNB       = cmdparams_get_int(&cmdparams,HCM_NB, NULL); //tuning position of NB Michelson correponding to center of Fe I line at rest (target wavelength)
00159   int HCMPOL      = cmdparams_get_int(&cmdparams,HCM_POL, NULL);//tuning position of tuning polarizer
00160   int N           = cmdparams_get_int(&cmdparams,NUM, NULL);    //number of tuning positions/wavelengths: 5 or 6
00161   int camera      = cmdparams_get_int(&cmdparams,khcamid, NULL);//front (1) or side (0) camera?
00162   int calibration = cmdparams_get_int(&cmdparams,"cal", NULL); 
00163 
00164   char COMMENT[256];
00165   strcpy(COMMENT,"Code used: lookup.c; CALIBRATION USED IS:"); //comment about the phase-map code
00166   if(calibration == 0)
00167     {
00168       printf("CALIBRATION USED IS CALIBRATION 11, VALID FROM MAY 2010 TO JANUARY 18, 2012\n");
00169       strcat(COMMENT," 11"); 
00170     }
00171   if(calibration == 1)
00172     {
00173       printf("CALIBRATION USED IS CALIBRATION 12, VALID FROM JANUARY 18, 2012\n");
00174       strcat(COMMENT," 12"); 
00175     }
00176   if(calibration == 2)
00177     {
00178       printf("CALIBRATION USED IS CALIBRATION 13, VALID FROM JANUARY 15, 2014\n");
00179       strcat(COMMENT," 13"); 
00180     }
00181   if(calibration != 0 && calibration != 1 && calibration != 2)
00182     {
00183       printf("the calibration requested does not exist\n");
00184       exit(1);
00185     }
00186   printf("COMMENT: %s\n",COMMENT);
00187 
00188   char *COMMENTS= "COMMENT";
00189 
00190   if(camera != 0 && camera != 1)
00191     {
00192       printf("Error: hcamid must be 0 or 1\n");
00193       exit(EXIT_FAILURE);
00194     }
00195   if(camera == 0) camera = LIGHT_SIDE;
00196   else            camera = LIGHT_FRONT;
00197   printf("CAMERA = %d\n",camera);
00198 
00199   if( N != 5 && N != 6 && N != 8 && N != 10)
00200     {
00201       printf("N must be equal to 5, 6, 8, or 10\n");
00202       exit(EXIT_FAILURE);
00203     }
00204 
00205   //TO CORRECT THE NON-TUNABLE TRANSMISSION PROFILE
00206   //****************************************************************************************
00207   float shiftw = 0.0; //0.01; //10 mA
00208   //*****************************************************************************************
00209 
00210   //HMI filter parameters
00211   float phases[nx2*ny2*5];                               //tunable elements, read from the input phase maps {nelement,nx2,ny2} (nx2=number of columns; ny2=number of rows)
00212 
00213   //wavelength grid
00214   int nlam      = 32001;//38441;//16001;                                 //the wavelength range must be, at least, [-8,+8] A around the target wavelength. That captures 99.76% of the filters integral values. MUST BE AN ODD NUMBER BECAUSE WE WANT THE ZERO
00215   double dlam   = dvtest*dlamdv;                         //wavelength resolution: depends on the resolution we want on the input velocities
00216   double lam[nlam];
00217 
00218   double FWHM,minimum;
00219   double ydefault = 1.;                                   //default value for the intensity of the solar line OUTSIDE a small range around 6173.3433 A (SHOULD BE 1)
00220   double ydefault2 = 0.;                                  //default value for the transmittance of the blocking filter and front window outside the range
00221   int i,j,k,ii,jj,iii,row,column,element;                 //loop increment variables
00222   FILE *fp;
00223 
00224   for(i=0;i<nlam;++i)    lam[i] = ((double)i-((double)nlam-1.0)/2.0)*dlam; //wavelength grid we will use
00225 
00226 
00227   //READING THE REFERENCE FeI LINE
00228   //WE WANT TO MAKE SURE THAT THE MINIMUM OF THE LINE IS AT THE TARGET WAVELENGTH 6173.3433 A
00229   //BECAUSE THIS MINIMUM DEFINES THE ZERO OF THE PHASES OF THE TUNABLE ELEMENTS
00230   //-----------------------------------------------------------------------------------------------------------------------------------
00231 
00232   float templineref[referencenlam];  //reference solar line: solar line of R. Ulrich at disk center, zero velocity, interpolated on a fine grid, and centered
00233   float wavelengthref[referencenlam];//wavelength grid for the reference wavelength
00234   double wavelength2[referencenlam];
00235   double templine0[nlam];            //will be the version of templineref interpolated onto the lam grid
00236 
00237   /*printf("Solar line profile file used = %s\n",referencesolarline);
00238   fp = fopen(referencesolarline,"rb");
00239   fread(templineref,sizeof(float),referencenlam,fp);
00240   fclose(fp);*/
00241 
00242   fp = fopen(referencewavelength,"rb");
00243   fread(wavelengthref,sizeof(float),referencenlam,fp);
00244   fclose(fp);
00245 
00246 
00247   //CAREFUL, I CHANGED THE WAVELENGTH GRID !!!!
00248   for(k=0;k<referencenlam;++k) wavelengthref[k]=(float)k/((float)referencenlam-1.0)-0.5; //varies between -0.5 and +0.5 A around the target wavelength (AND THERE IS NO ZERO, WHICH IS USEFUL BECAUSE OTHERWISE THERE WOULD BE A NAN IN THE VOIGT PROFILE)
00249 
00250 
00251   /*
00252   //find the minimum of the line
00253   minimum=ydefault;
00254   k=0;
00255   for(j=0;j<referencenlam;j++) 
00256     {   
00257       if(templineref[j] < minimum)
00258     {
00259       minimum=templineref[j];
00260       k=j;
00261     }
00262     }
00263     printf("THE MINIMUM OF THE SOLAR LINE IS AT LAM= %f\n",wavelengthref[k]);*/
00264 
00265   //parameters of the Fe I line profiles observed by Roger Ulrich (OBTAINED BY CODE centertolimb.pro)
00266   //-------------------------------------------------------------------------------------------------
00267   //double cost[3]={1.0,0.70710678,0.5};                    //cos(theta) where theta is the angular distance from disk center
00268   double minimumCoeffs[2]={0.41922611,0.24190794};        //minimum intensity (Id if I=1-Id*exp() ), assuming the continuum is at 1 (result from a Gaussian fit in the range [-0.055,0.055] Angstroms)
00269   double FWHMCoeffs[2]={151.34559,-58.521771};            //FWHM of the solar line in milliAngstrom (result from a Gaussian fit in the range [-0.055,0.055] Angstroms)
00270 
00271   //ALTERNATIVE Fe I LINES
00272   //--------------------------------------------------------------------------------------------------------------------------------
00273   //printf("!!!!!!!! WARNING !!!!!!!!\n");
00274   //printf("WE ARE USING A GAUSSIAN LINE PROFILE !\n");
00275   float Icg;
00276   float fdepthg;
00277   float fwidthg;
00278 
00279   if(calibration == 0)
00280     {
00281       Icg=1.0;
00282       fdepthg=0.5625;//0.63;//0.56;//0.642710;//0.651742;
00283       fwidthg=0.06415;//0.069334;//0.071124;//0.0687;//0.069334;//0.069598;
00284     }
00285   if(calibration == 1)
00286     {
00287       Icg=1.0;
00288       fdepthg=0.53;
00289       fwidthg=0.0615;
00290     }
00291   if(calibration == 2)
00292     {
00293       Icg=1.0;
00294       fdepthg=0.58;
00295       fwidthg=0.058;
00296     }
00297 
00298   //1 Gaussian
00299   //for(k=0;k<referencenlam;++k) templineref[k]=Icg-fdepthg*exp(-wavelengthref[k]*wavelengthref[k]/fwidthg/fwidthg);
00300 
00301 
00302   //2 Gaussians (1 big and 1 small to simulate the asymmetry)
00303   //for(k=0;k<referencenlam;++k) templineref[k]=Icg-fdepthg*exp(-wavelengthref[k]*wavelengthref[k]/fwidthg/fwidthg)-0.025663363*exp(-(wavelengthref[k]+2.*0.0463867)*(wavelengthref[k]+2.*0.0463867)/0.18869677/0.18869677);
00304 
00305 
00306   //Voigt profile
00307   float a=0.03; //damping
00308   if(calibration == 2) a=-0.09;
00309   float l;
00310 
00311   /*for(k=0;k<referencenlam;++k) 
00312     {
00313       l=wavelengthref[k]/fwidthg;
00314       if(fabs(wavelengthref[k]) <= 1.5) templineref[k]= Icg-fdepthg*exp(-l*l)*(1.0-a*2.0/sqrt(M_PI)*(1.0/2.0/l/l)*((4.0*l*l+3.0)*(l*l+1.0)*exp(-l*l)-1.0/l/l*(2.0*l*l+3.0)*sinh(l*l)));
00315       else templineref[k]=Icg;
00316 
00317       }*/
00318 
00319 
00320   //Voigt profile + Gaussian profile for the small Lanthanum II line at -0.615 A from target wavelength
00321   /*for(k=0;k<referencenlam;++k) 
00322     {
00323       l=wavelengthref[k]/fwidthg;
00324       if(fabs(wavelengthref[k]) <= 1.5) templineref[k]= Icg-fdepthg*exp(-l*l)*(1.0-a*2.0/sqrt(M_PI)*(1.0/2.0/l/l)*((4.0*l*l+3.0)*(l*l+1.0)*exp(-l*l)-1.0/l/l*(2.0*l*l+3.0)*sinh(l*l)))-0.0225*exp(-(wavelengthref[k]+0.615)*(wavelengthref[k]+0.615)/0.076/0.076);
00325       else templineref[k]=Icg;
00326       
00327       }*/
00328 
00329   //Voigt profile + Gaussian profiles in the wing to simulate asymmetry
00330 
00331   if(calibration == 0)
00332     {
00333       for(k=0;k<referencenlam;++k) 
00334     {
00335       l=wavelengthref[k]/fwidthg;
00336       if(fabs(l) <= 26.5) templineref[k]= Icg-fdepthg*exp(-l*l)*(1.0-a*2.0/sqrt(M_PI)*(1.0/2.0/l/l)*((4.0*l*l+3.0)*(l*l+1.0)*exp(-l*l)-1.0/l/l*(2.0*l*l+3.0)*sinh(l*l)))-0.015*exp(-(wavelengthref[k]+0.225)*(wavelengthref[k]+0.225)/0.2/0.2)+0.004*exp(-(wavelengthref[k]-0.150)*(wavelengthref[k]-0.150)/0.22/0.22); //sinh((x>26.5)^2) = NAN
00337       else templineref[k]=Icg-0.015*exp(-(wavelengthref[k]+0.225)*(wavelengthref[k]+0.225)/0.2/0.2)+0.004*exp(-(wavelengthref[k]-0.150)*(wavelengthref[k]-0.150)/0.22/0.22);
00338     }
00339     }
00340   if(calibration == 1)
00341     {
00342       for(k=0;k<referencenlam;++k) 
00343     {
00344       l=wavelengthref[k]/fwidthg;
00345       if(fabs(l) <= 26.5) templineref[k]= Icg-fdepthg*exp(-l*l)*(1.0-a*2.0/sqrt(M_PI)*(1.0/2.0/l/l)*((4.0*l*l+3.0)*(l*l+1.0)*exp(-l*l)-1.0/l/l*(2.0*l*l+3.0)*sinh(l*l)))+0.01*exp(-(wavelengthref[k]+0.225)*(wavelengthref[k]+0.225)/0.2/0.2)+0.015*exp(-(wavelengthref[k]-0.1)*(wavelengthref[k]-0.1)/0.25/0.25); //sinh((x>26.5)^2) = NAN
00346       else templineref[k]=Icg+0.01*exp(-(wavelengthref[k]+0.225)*(wavelengthref[k]+0.225)/0.2/0.2)+0.015*exp(-(wavelengthref[k]-0.1)*(wavelengthref[k]-0.1)/0.25/0.25);
00347     }
00348     }
00349   if(calibration == 2)
00350     {
00351       for(k=0;k<referencenlam;++k) 
00352     {
00353       l=wavelengthref[k]/fwidthg;
00354       if(fabs(l) <= 26.5) templineref[k]= Icg-fdepthg*exp(-l*l)*(1.0-a*2.0/sqrt(M_PI)*(1.0/2.0/l/l)*((4.0*l*l+3.0)*(l*l+1.0)*exp(-l*l)-1.0/l/l*(2.0*l*l+3.0)*sinh(l*l)))+0.0074*exp(-(wavelengthref[k]+0.2)*(wavelengthref[k]+0.2)/0.13/0.13)+0.021*exp(-(wavelengthref[k]-0.05)*(wavelengthref[k]-0.05)/0.18/0.18); //sinh((x>26.5)^2) = NAN
00355       else templineref[k]=Icg+0.0074*exp(-(wavelengthref[k]+0.2)*(wavelengthref[k]+0.2)/0.13/0.13)+0.021*exp(-(wavelengthref[k]-0.05)*(wavelengthref[k]-0.05)/0.18/0.18);
00356     }
00357     }
00358 
00359 
00360 
00361 
00362   //for(k=0;k<referencenlam;++k) templineref[k]=1.0-phases[4+location*5]*exp(-wavelengthref[k]*wavelengthref[k]/phases[3+location*5]/phases[3+location*5]);
00363   //find the minimum of the line
00364   minimum=ydefault;
00365   k=0;
00366   for(j=0;j<referencenlam;j++) 
00367     {   
00368       if(templineref[j] < minimum)
00369     {
00370       minimum=templineref[j];
00371       k=j;
00372     }
00373     }
00374   printf("THE MINIMUM OF THE SOLAR LINE IS AT LAM= %f\n",wavelengthref[k]);
00375 
00376   //parameters for a Gaussian profile
00377   //-------------------------------------------------------------------------------------------------
00378   //double minimumCoeffs=0.515347;                          //linedepth, not minimum intensity, at disk center
00379   //double FWHMCoeffs=100.67102;                            //FWHM at disk center
00380 
00381   const char *HCMNBs      = "HWL4POS";                   //keyword for the position of the HCM of NB Michelson
00382   const char *HCMPOLs     = "HWL3POS";                   //keyword for the position of the HCM of the tuning polarizer 
00383   const char *HCMWBs      = "HWL2POS";                   //keyword for the position of the HCM of WB Michelson
00384   const char *HCME1s      = "HWL1POS";                   //keyword for the position of the HCM of Lyot E1
00385   const char *Ns          = "NWL";                       //keyword for the number of tuning positions
00386   const char *keyname     = "FSN_START";     
00387   const char *keyname2    = "T_START";       
00388   const char *keyname3    = "FSN_STOP";      
00389   const char *keyname4    = "T_STOP";      
00390   const char *keyname5    = "T_REC";
00391   const char *keyname6    = "FSN_REC";
00392   const char *keyname7    = "HCAMID";
00393   const char *RADIUS      = "R_SUN";                      //maximum radius for which the phases of the tunable elements where computed (in pixels)
00394   const char *X0c         = "CRPIX1";                     //x-axis location of image disk in pixels, starting at 1
00395   const char *Y0c         = "CRPIX2";                     //y-axis location of image disk in pixels
00396   const char *SCALE       = "CDELT1";                     //image scale in the x direction
00397   const char *CROTA2      = "CROTA2";                     // added for ROLL 2018
00398   const char *OBS_VR      = "OBS_VR";                     // added for ROLL 2018
00399   const char *HCMNB0      = "HCMNB";
00400   const char *HCMWB0      = "HCMWB";
00401   const char *HCME10      = "HCME1";
00402   const char *FSRNBs      = "FSRNB";
00403   const char *FSRWBs      = "FSRWB";
00404   const char *FSRE1s      = "FSRE1";
00405   const char *FSRE2s      = "FSRE2";
00406   const char *FSRE3s      = "FSRE3";
00407   const char *FSRE4s      = "FSRE4";
00408   const char *FSRE5s      = "FSRE5";
00409   const char *CBLOCKERs   = "CBLOCKER";
00410   const char *DATEs       = "DATE";
00411 
00412 
00413   double inttune = (N-1.)/2.;                             //number of tuning intervals on each wing of the solar line
00414   double dtune   = FSR[0]/2.5;                            //nominal wavelength separation between each tuning position: we use 2.5 instead of inttune so that the spacing is the same for 5 or 6 positions
00415   double dvtune  = dtune/dlamdv;
00416   double *angle  = NULL;
00417 
00418   double X0 = (nx2-1)/2.0;                                //center of the "solar disk" in pixel units (starts at 0.0)
00419   double Y0 = (ny2-1)/2.0;
00420   int  ratio= 4096/nx2;
00421   float solarradiusmax = 2047.5-(float)ratio;             //maximum radius of computation of the look-up tables (in pixels), we subtract ratio pixels to avoid edge issues
00422 
00423   angle=(double *)malloc(N*sizeof(double));
00424   if(angle == NULL)
00425     {
00426       printf("Error: unable to allocate memory to angle\n");
00427       exit(EXIT_FAILURE);
00428     }
00429 
00430   printf("Number of tuning positions= %d\n",N);
00431 
00432   if(N == 6)
00433     {
00434       angle[0]=+2.5;//tuning positions in the order I0, I1, I2, I3, I4, and I5 (FROM RED WING TO BLUE WING)
00435       angle[1]=+1.5;
00436       angle[2]=+0.5;
00437       angle[3]=-0.5;
00438       angle[4]=-1.5;
00439       angle[5]=-2.5;
00440     }
00441   if(N == 5)
00442     {
00443       angle[0]=+2.0;
00444       angle[1]=+1.0;
00445       angle[2]= 0.0;
00446       angle[3]=-1.0;
00447       angle[4]=-2.0;
00448     }
00449   if(N == 8)
00450     {
00451       angle[7]=+3.5;//tuning positions in the order I7, I0, I1, I2, I3, I4, I5, and I6 (RANDOM ORDER !)
00452       angle[0]=+2.5;
00453       angle[1]=+1.5;
00454       angle[2]=+0.5;
00455       angle[3]=-0.5;
00456       angle[4]=-1.5;
00457       angle[5]=-2.5;
00458       angle[6]=-3.5;
00459       ntest=1333 ;//WE CHANGE THE NUMBER OF TEST VELOCITIES BECAUSE THE DYNAMICAL RANGE IS INCREASED (OVERRIDE WHAT IS IN HMIparam.h). MUST BE AN ODD NUMBER
00460     }
00461   if(N == 10)
00462     {
00463       angle[9]=+4.5;//tuning positions in the order I9, I7, I0, I1, I2, I3, I4, I5, I6, and I8 (RANDOM ORDER !)
00464       angle[7]=+3.5;
00465       angle[0]=+2.5;
00466       angle[1]=+1.5;
00467       angle[2]=+0.5;
00468       angle[3]=-0.5;
00469       angle[4]=-1.5;
00470       angle[5]=-2.5;
00471       angle[6]=-3.5;
00472       angle[8]=-4.5;
00473       ntest=1333 ;//WE CHANGE THE NUMBER OF TEST VELOCITIES BECAUSE THE DYNAMICAL RANGE IS INCREASED (OVERRIDE WHAT IS IN HMIparam.h). MUST BE AN ODD NUMBER
00474     }
00475 
00476 
00477   double pv1 = 0.0;
00478   double pv2 = 0.0;
00479 
00480   double lineint[nlam];
00481   double f1c = 0.; //First Fourier coefficients
00482   double f1s = 0.;
00483   double f2c = 0.; //Second Fourier coefficients
00484   double f2s = 0.;
00485   //double norm[N];  //norm of the filters
00486   double distance[nx2*ny2];
00487 
00488   double WRONGDISTANCE=100.;
00489   double BUFFERDISTANCE=0.0; //value of cos(theta) in the buffer zone
00490 
00491   double *cosi = NULL;
00492   double *sini = NULL;
00493   double *cos2i= NULL;
00494   double *sin2i= NULL;
00495   cosi = (double *)malloc(N*sizeof(double));  
00496   sini = (double *)malloc(N*sizeof(double));  
00497   cos2i= (double *)malloc(N*sizeof(double));  
00498   sin2i= (double *)malloc(N*sizeof(double));  
00499   if(cosi == NULL || sini == NULL || cos2i == NULL || sin2i == NULL)
00500     {
00501       printf("Error: memory could not be allocated to the Fourier coefficient arrays\n");
00502       exit(EXIT_FAILURE);
00503     }
00504   for(i=0;i<N;++i) 
00505     {
00506       angle[i]= angle[i]*2.0*M_PI/(double)N;
00507       cosi[i] = cos(angle[i]);
00508       sini[i] = sin(angle[i]);
00509       cos2i[i]= cos(2.0*angle[i]);
00510       sin2i[i]= sin(2.0*angle[i]);
00511     }
00512   
00513   //printf("NTEST= %d\n",ntest);
00514 
00515   int axis[3]    = {5,nx2,ny2};                                   //size of input phase maps of tunable elements (nx2= number of columns, ny2= number of rows)
00516   int axisout[3] = {ntest*2,nx2,ny2};                             //size of the output look-up tables
00517 
00518   if((axisout[0] != ntest*2) || (axisout[1] > maxNx) || (axisout[2] > maxNx) ||(axisout[1] != axisout[2]) )
00519     {
00520       printf("Error: dimensions are incorrect %d %d %d\n",axisout[0],axisout[1],axisout[2]);
00521       exit(EXIT_FAILURE);
00522     }
00523 
00524   int  status = DRMS_SUCCESS; 
00525   DRMS_Array_t *arrout   = NULL;                                  //array that will contain the phase maps 
00526   DRMS_Type_t type;
00527   type    = DRMS_TYPE_FLOAT;                                      //type of the output data: FLOAT (NOT DOUBLE, TO SAVE I/O TIME WITH OBSERVABLE CODE)
00528   arrout  = drms_array_create(type,3,axisout,NULL,&status);       //create the array that will contain the phase maps
00529   if(status != DRMS_SUCCESS)
00530     {
00531       printf("Error: unable to create array arrout\n");
00532       exit(EXIT_FAILURE);
00533     }
00534   float *vel  = arrout->data;                                     //lookup tables (VELOCITIES RETURNED BY MDI-LIKE ALGORITHM AS A FUNCTION OF INPUT VELOCITIES)
00535   memset(arrout->data,0.0,drms_array_size(arrout));               //initializes the look-up table values to 0.0
00536 
00537   double vtest[ntest];
00538 
00539   for(i=0;i<ntest;++i) vtest[i] = dvtest*((double)i-((double)ntest-1.0)/2.0);
00540   printf("VTEST MIN AND MAX = %f %f\n",vtest[0],vtest[ntest-1]);
00541   pv1 = dvtune*(double)(N-1); //WARNING THIS IS WRONG, SHOULD BE dvtune*(double)N
00542   pv2 = pv1/2.;
00543 
00544   /***********************************************************************************************************/
00545   /*INPUT SERIES                                                                                           */
00546   /***********************************************************************************************************/
00547 
00548   char *inRecQuery  = cmdparams_get_str(&cmdparams, kRecSetIn, NULL); //cmdparams is defined in jsoc_main.h
00549   DRMS_RecordSet_t *data = drms_open_records(drms_env,inRecQuery,&status);   //open the records from the input series
00550   int nRecs;
00551 
00552   if (status == DRMS_SUCCESS && data != NULL && data->n > 0)
00553     {
00554       nRecs = data->n;                                           //number of records in the input series 
00555       printf("Number of phase-map satisfying the request= %d \n",nRecs);    //check if the number of records is appropriate
00556     }
00557   else
00558     {
00559       printf("Input phase-map series %s doesn't exist\n",inRecQuery);//if the input series does not exit
00560       exit(EXIT_FAILURE);                                            //we exit the program
00561     }
00562 
00563       DRMS_Segment_t *segin  = NULL;
00564       DRMS_Array_t *arrin    = NULL;
00565       int  keyvalue          = 0;
00566       DRMS_Type_Value_t  keyvalue2;
00567       int  keyvalue3         = 0;
00568       DRMS_Type_Value_t  keyvalue4;
00569       int keyvalue5          = 0;
00570       DRMS_Type_Value_t  keyvalue6;
00571       int keyvalue7;
00572       int recofinterest;
00573       DRMS_Type_t type2      = DRMS_TYPE_TIME;
00574       float *tempphase = NULL; 
00575       TIME   interntime;
00576       TIME   interntime2;
00577       TIME   interntime3;
00578       float XCENTER,YCENTER,RSUN,CDELT1;
00579       double crota2, obs_vr;            // added for ROLL 2018
00580 
00581       //LOCATE THE PHASE-MAP RECORD TO USE
00582       i=0;
00583       while(i<nRecs)
00584     {
00585       keyvalue7   = drms_getkey_int(data->records[i],keyname7,&status);
00586       if(status != DRMS_SUCCESS)
00587         {
00588           printf("Error: unable to read keyword %s\n",keyname7);
00589           exit(EXIT_FAILURE);
00590         }
00591       if(keyvalue7 == camera)
00592         {
00593           recofinterest=i; //the record contains the phase maps obtained from data taken with the camera we want
00594           break;
00595         }
00596       i++;
00597     }
00598       if(i == nRecs)
00599     {
00600       printf("Error: no phase-map record was found with the proper FSN_REC and HCAMID\n");
00601       exit(EXIT_FAILURE);
00602     }
00603 
00604       printf("RECORD OF INTEREST = %d\n",i);
00605 
00606       segin     = drms_segment_lookupnum(data->records[recofinterest],0);
00607       arrin     = drms_segment_read(segin,segin->info->type,&status);
00608       if(status != DRMS_SUCCESS)
00609     {
00610       printf("Error: unable to read a data segment\n");
00611       exit(EXIT_FAILURE);
00612     }
00613       printf("Segment read\n");
00614       tempphase = arrin->data;
00615       //for(i=0;i<nx2*ny2*nelement;++i) phases[i] = tempphase[i]*M_PI/180.; //NB: PHASE MAPS ARE ASSUMED TO BE IN TYPE FLOAT AND IN DEGREES
00616       for(i=0;i<nx2*ny2*5;++i) phases[i] = tempphase[i]; //NB: PHASE MAPS ARE ASSUMED TO BE IN TYPE FLOAT AND IN DEGREES
00617 
00618 
00619       printf("Reading keywords\n");
00620       //READ A FEW KEYWORDS FROM THE PHASE-MAP SERIES
00621       XCENTER         = drms_getkey_float(data->records[recofinterest],X0c,&status);
00622       if(status != DRMS_SUCCESS)
00623     {
00624       printf("Error: unable to read the keyword %s\n",X0c);
00625       exit(EXIT_FAILURE);
00626     }
00627       else XCENTER=XCENTER-1.0; //BECAUSE CRPIX1 STARTS AT 1 
00628       YCENTER         = drms_getkey_float(data->records[recofinterest],Y0c,&status);
00629       if(status != DRMS_SUCCESS)
00630     {
00631       printf("Error: unable to read the keyword %s\n",Y0c);
00632       exit(EXIT_FAILURE);
00633     }
00634       else YCENTER=YCENTER-1.0; //BECAUSE CRPIX2 STARTS AT 1
00635       CDELT1  = drms_getkey_float(data->records[recofinterest],SCALE,&status);
00636       if(status != DRMS_SUCCESS)
00637     {
00638       printf("Error: unable to read the keyword %s\n",SCALE);
00639       exit(EXIT_FAILURE);
00640     }
00641       RSUN  = drms_getkey_float(data->records[recofinterest],RADIUS,&status);
00642       if(status != DRMS_SUCCESS)
00643     {
00644       printf("Error: unable to read the keyword %s\n",RADIUS);
00645       exit(EXIT_FAILURE);
00646     }
00647 
00648       //CHECK THAT THESE KEYWORDS ARE WITHIN ACCEPTABLE LIMITS
00649       if(fabs(XCENTER-2047.5) > ratio/2)
00650     {
00651       printf("Error: the value of %s of the phase-map series is too far from the center of the CCD\n",X0c);
00652       exit(EXIT_FAILURE);
00653     }
00654       if(fabs(YCENTER-2047.5) > ratio/2)
00655     {
00656       printf("Error: the value of %s of the phase-map series is too far from the center of the CCD\n",Y0c);
00657       exit(EXIT_FAILURE);
00658     }
00659       if(fabs(RSUN-solarradiusmax) > ratio/2)
00660     {
00661       printf("Error: the radius %s of the phase-map series is too different from %f\n",RADIUS,solarradiusmax);
00662       exit(EXIT_FAILURE);
00663     } 
00664       printf("Keywords read\n");
00665 
00666 
00667       //CALCULATE DISTANCES
00668       for(i=0;i<nx2*ny2;++i) 
00669     {
00670       row    = i / nx2;//nx2= number of columns
00671       column = i % nx2;
00672       distance[i]=sqrt(((double)column-X0)*((double)column-X0)+((double)row-Y0)*((double)row-Y0))*NOMINALSCALE*(double)ratio; //radial distance from disk center in arcseconds
00673       if(distance[i] <= (double)solarradiustable) distance[i]=cos(asin(distance[i]/(double)solarradiustable));                //convert projected angular distance from disk center to cos(angle) between the solar surface normal and the l.o.s. to the observer (INFINITE DISTANCE APPROXIMATION)
00674       else
00675         {
00676           if(distance[i] <= (double)solarradiusmax*NOMINALSCALE) distance[i]=BUFFERDISTANCE;//to create a buffer zone
00677           else distance[i]=WRONGDISTANCE;
00678         }
00679     }
00680 
00681       keyvalue   = drms_getkey_int(data->records[recofinterest],keyname,&status);
00682       if(status != DRMS_SUCCESS)
00683     {
00684       printf("Error: unable to read keyword %s\n",keyname);
00685       exit(EXIT_FAILURE);
00686     }
00687       keyvalue2  = drms_getkey(data->records[recofinterest],keyname2,&type2,&status);
00688       if(status != DRMS_SUCCESS)
00689     {
00690       printf("Error: unable to read keyword %s\n",keyname2);
00691       exit(EXIT_FAILURE);
00692     }
00693       keyvalue3  = drms_getkey_int(data->records[recofinterest],keyname3,&status);
00694       if(status != DRMS_SUCCESS)
00695     {
00696       printf("Error: unable to read keyword %s\n",keyname3);
00697       exit(EXIT_FAILURE);
00698     }
00699       keyvalue4  = drms_getkey(data->records[recofinterest],keyname4,&type2,&status);
00700       if(status != DRMS_SUCCESS)
00701     {
00702       printf("Error: unable to read keyword %s\n",keyname4);
00703       exit(EXIT_FAILURE);
00704     }
00705       keyvalue5  = drms_getkey_int(data->records[recofinterest],keyname6,&status);
00706       if(status != DRMS_SUCCESS)
00707     {
00708       printf("Error: unable to read keyword %s\n",keyname6);
00709       exit(EXIT_FAILURE);
00710     }
00711       keyvalue6  = drms_getkey(data->records[recofinterest],keyname5,&type2,&status);
00712       if(status != DRMS_SUCCESS)
00713     {
00714       printf("Error: unable to read keyword %s\n",keyname5);
00715       exit(EXIT_FAILURE);
00716     }
00717       printf("FSN_START value %d\n",keyvalue);                         
00718       interntime = keyvalue2.time_val;
00719       interntime2= keyvalue4.time_val;
00720       interntime3= keyvalue6.time_val;
00721 
00722 
00723   /***********************************************************************************************************/
00724   /*INFORMATION ON THE COTUNE SEQUENCE USED                                                                  */
00725   /***********************************************************************************************************/
00726 
00727       //the user provides the tuning position in HCM steps for the central tuning position of the framelist (i.e.
00728       //the tuning positions for which the phases of the elements are closest to their average value over the
00729       //phase maps
00730 
00731       double *HCME1phase,*HCMWBphase,*HCMNBphase;
00732       HCME1phase = (double *) malloc(N*sizeof(double));
00733       HCMWBphase = (double *) malloc(N*sizeof(double));
00734       HCMNBphase = (double *) malloc(N*sizeof(double));
00735 
00736 
00737   /***********************************************************************************************************/
00738   /*TUNING POSITIONS                                                                                         */
00739   /***********************************************************************************************************/
00740 
00741 
00742      if(N == 6)
00743        {
00744      HCME1phase[0]= (double) ((HCME1+15)*6 % 360)*M_PI/180.0; //I0
00745      HCME1phase[1]= (double) ((HCME1+9 )*6 % 360)*M_PI/180.0; //I1
00746      HCME1phase[2]= (double) ((HCME1+3 )*6 % 360)*M_PI/180.0; //I2
00747      HCME1phase[3]= (double) ((HCME1-3 )*6 % 360)*M_PI/180.0; //I3
00748      HCME1phase[4]= (double) ((HCME1-9 )*6 % 360)*M_PI/180.0; //I4
00749      HCME1phase[5]= (double) ((HCME1-15)*6 % 360)*M_PI/180.0; //I5
00750                           
00751      HCMWBphase[0]= (double) ((HCMWB-30)*6 % 360)*M_PI/180.0;
00752      HCMWBphase[1]= (double) ((HCMWB-18)*6 % 360)*M_PI/180.0;
00753      HCMWBphase[2]= (double) ((HCMWB-6 )*6 % 360)*M_PI/180.0;
00754      HCMWBphase[3]= (double) ((HCMWB+6 )*6 % 360)*M_PI/180.0;
00755      HCMWBphase[4]= (double) ((HCMWB+18)*6 % 360)*M_PI/180.0;
00756      HCMWBphase[5]= (double) ((HCMWB-30)*6 % 360)*M_PI/180.0;
00757                           
00758      HCMNBphase[0]= (double) ((HCMNB-0 )*6 % 360)*M_PI/180.0;
00759      HCMNBphase[1]= (double) ((HCMNB+24)*6 % 360)*M_PI/180.0;
00760      HCMNBphase[2]= (double) ((HCMNB-12)*6 % 360)*M_PI/180.0;
00761      HCMNBphase[3]= (double) ((HCMNB+12)*6 % 360)*M_PI/180.0;
00762      HCMNBphase[4]= (double) ((HCMNB-24)*6 % 360)*M_PI/180.0;
00763      HCMNBphase[5]= (double) ((HCMNB+0 )*6 % 360)*M_PI/180.0;
00764        }
00765 
00766      if(N == 5)
00767        {
00768 
00769      HCME1phase[0]= (double) ((HCME1+12)*6 % 360)*M_PI/180.0; //I0
00770      HCME1phase[1]= (double) ((HCME1+6 )*6 % 360)*M_PI/180.0; //I1
00771      HCME1phase[2]= (double) ((HCME1+0 )*6 % 360)*M_PI/180.0; //I2
00772      HCME1phase[3]= (double) ((HCME1-6 )*6 % 360)*M_PI/180.0; //I3
00773      HCME1phase[4]= (double) ((HCME1-12)*6 % 360)*M_PI/180.0; //I4
00774                           
00775      HCMWBphase[0]= (double) ((HCMWB-24)*6 % 360)*M_PI/180.0;
00776      HCMWBphase[1]= (double) ((HCMWB-12)*6 % 360)*M_PI/180.0;
00777      HCMWBphase[2]= (double) ((HCMWB+0 )*6 % 360)*M_PI/180.0;
00778      HCMWBphase[3]= (double) ((HCMWB+12)*6 % 360)*M_PI/180.0;
00779      HCMWBphase[4]= (double) ((HCMWB+24)*6 % 360)*M_PI/180.0;
00780                           
00781      HCMNBphase[0]= (double) ((HCMNB+12)*6 % 360)*M_PI/180.0;
00782      HCMNBphase[1]= (double) ((HCMNB-24)*6 % 360)*M_PI/180.0;
00783      HCMNBphase[2]= (double) ((HCMNB+0 )*6 % 360)*M_PI/180.0;
00784      HCMNBphase[3]= (double) ((HCMNB+24)*6 % 360)*M_PI/180.0;
00785      HCMNBphase[4]= (double) ((HCMNB-12)*6 % 360)*M_PI/180.0;
00786        }
00787 
00788 
00789      if(N == 8)
00790        {
00791      HCME1phase[7]= (double) ((HCME1+21)*6 % 360)*M_PI/180.0; //I7
00792      HCME1phase[0]= (double) ((HCME1+15)*6 % 360)*M_PI/180.0; //I0
00793      HCME1phase[1]= (double) ((HCME1+9 )*6 % 360)*M_PI/180.0; //I1
00794      HCME1phase[2]= (double) ((HCME1+3 )*6 % 360)*M_PI/180.0; //I2
00795      HCME1phase[3]= (double) ((HCME1-3 )*6 % 360)*M_PI/180.0; //I3
00796      HCME1phase[4]= (double) ((HCME1-9 )*6 % 360)*M_PI/180.0; //I4
00797      HCME1phase[5]= (double) ((HCME1-15)*6 % 360)*M_PI/180.0; //I5
00798      HCME1phase[6]= (double) ((HCME1-21)*6 % 360)*M_PI/180.0; //I6
00799                         
00800      HCMWBphase[7]= (double) ((HCMWB+18)*6 % 360)*M_PI/180.0;
00801      HCMWBphase[0]= (double) ((HCMWB-30)*6 % 360)*M_PI/180.0;
00802      HCMWBphase[1]= (double) ((HCMWB-18)*6 % 360)*M_PI/180.0;
00803      HCMWBphase[2]= (double) ((HCMWB-6 )*6 % 360)*M_PI/180.0;
00804      HCMWBphase[3]= (double) ((HCMWB+6 )*6 % 360)*M_PI/180.0;
00805      HCMWBphase[4]= (double) ((HCMWB+18)*6 % 360)*M_PI/180.0;
00806      HCMWBphase[5]= (double) ((HCMWB-30)*6 % 360)*M_PI/180.0;
00807      HCMWBphase[6]= (double) ((HCMWB-18)*6 % 360)*M_PI/180.0;
00808 
00809      HCMNBphase[7]= (double) ((HCMNB-24)*6 % 360)*M_PI/180.0;
00810      HCMNBphase[0]= (double) ((HCMNB-0 )*6 % 360)*M_PI/180.0;
00811      HCMNBphase[1]= (double) ((HCMNB+24)*6 % 360)*M_PI/180.0;
00812      HCMNBphase[2]= (double) ((HCMNB-12)*6 % 360)*M_PI/180.0;
00813      HCMNBphase[3]= (double) ((HCMNB+12)*6 % 360)*M_PI/180.0;
00814      HCMNBphase[4]= (double) ((HCMNB-24)*6 % 360)*M_PI/180.0;
00815      HCMNBphase[5]= (double) ((HCMNB+0 )*6 % 360)*M_PI/180.0;
00816      HCMNBphase[6]= (double) ((HCMNB+24 )*6 % 360)*M_PI/180.0;
00817        }
00818 
00819      if(N == 10)
00820        {
00821      HCME1phase[9]= (double) ((HCME1+27)*6 % 360)*M_PI/180.0; //I9
00822      HCME1phase[7]= (double) ((HCME1+21)*6 % 360)*M_PI/180.0; //I7
00823      HCME1phase[0]= (double) ((HCME1+15)*6 % 360)*M_PI/180.0; //I0
00824      HCME1phase[1]= (double) ((HCME1+9 )*6 % 360)*M_PI/180.0; //I1
00825      HCME1phase[2]= (double) ((HCME1+3 )*6 % 360)*M_PI/180.0; //I2
00826      HCME1phase[3]= (double) ((HCME1-3 )*6 % 360)*M_PI/180.0; //I3
00827      HCME1phase[4]= (double) ((HCME1-9 )*6 % 360)*M_PI/180.0; //I4
00828      HCME1phase[5]= (double) ((HCME1-15)*6 % 360)*M_PI/180.0; //I5
00829      HCME1phase[6]= (double) ((HCME1-21)*6 % 360)*M_PI/180.0; //I6
00830      HCME1phase[8]= (double) ((HCME1-27)*6 % 360)*M_PI/180.0; //I8
00831                         
00832      HCMWBphase[9]= (double) ((HCMWB+6)*6 % 360)*M_PI/180.0;
00833      HCMWBphase[7]= (double) ((HCMWB+18)*6 % 360)*M_PI/180.0;
00834      HCMWBphase[0]= (double) ((HCMWB-30)*6 % 360)*M_PI/180.0;
00835      HCMWBphase[1]= (double) ((HCMWB-18)*6 % 360)*M_PI/180.0;
00836      HCMWBphase[2]= (double) ((HCMWB-6 )*6 % 360)*M_PI/180.0;
00837      HCMWBphase[3]= (double) ((HCMWB+6 )*6 % 360)*M_PI/180.0;
00838      HCMWBphase[4]= (double) ((HCMWB+18)*6 % 360)*M_PI/180.0;
00839      HCMWBphase[5]= (double) ((HCMWB-30)*6 % 360)*M_PI/180.0;
00840      HCMWBphase[6]= (double) ((HCMWB-18)*6 % 360)*M_PI/180.0;
00841      HCMWBphase[8]= (double) ((HCMWB-6)*6 % 360)*M_PI/180.0;
00842 
00843      HCMNBphase[9]= (double) ((HCMNB+12 )*6 % 360)*M_PI/180.0;
00844      HCMNBphase[7]= (double) ((HCMNB-24)*6 % 360)*M_PI/180.0;
00845      HCMNBphase[0]= (double) ((HCMNB-0 )*6 % 360)*M_PI/180.0;
00846      HCMNBphase[1]= (double) ((HCMNB+24)*6 % 360)*M_PI/180.0;
00847      HCMNBphase[2]= (double) ((HCMNB-12)*6 % 360)*M_PI/180.0;
00848      HCMNBphase[3]= (double) ((HCMNB+12)*6 % 360)*M_PI/180.0;
00849      HCMNBphase[4]= (double) ((HCMNB-24)*6 % 360)*M_PI/180.0;
00850      HCMNBphase[5]= (double) ((HCMNB+0 )*6 % 360)*M_PI/180.0;
00851      HCMNBphase[6]= (double) ((HCMNB+24 )*6 % 360)*M_PI/180.0;
00852      HCMNBphase[8]= (double) ((HCMNB-12)*6 % 360)*M_PI/180.0;
00853        }
00854 
00855 
00856 
00857   /***********************************************************************************************************/
00858   /*OUTPUT SERIES                                                                                            */
00859   /***********************************************************************************************************/
00860 
00861   char *dsout = cmdparams_get_str(&cmdparams, kDSOut, NULL);         //series name of the output look-up table series
00862   drms_series_exists(drms_env, dsout, &status);                      //check whether or not output series exists
00863                                                                      //drms_env is defined in jsoc_main.h
00864   if (status == DRMS_ERROR_UNKNOWNSERIES)
00865     {
00866       printf("Output Lookup Table series %s doesn't exist\n",dsout);  //if the output series does not exit
00867       exit(EXIT_FAILURE);                                            //we exit the program
00868     } 
00869   if (status == DRMS_SUCCESS)
00870     {
00871       printf("Output Lookup Table series %s exists.\n",dsout);
00872     }
00873 
00874   /***********************************************************************************************************/
00875   /*BLOCKER FILTER + FRONT WINDOW                                                                            */
00876   /***********************************************************************************************************/
00877 
00878   printf("Reading the front window and blocker filter spatially-averaged transmission profiles\n");
00879   for(i=0;i<nfront;++i)
00880     {
00881       wavelengthd[i]=wavelengthd[i]*10.0-lam0;
00882       frontwindowd[i]=frontwindowd[i]/100.0;
00883     }
00884   for(i=0;i<nblocker;++i)
00885     {
00886       wavelengthbd[i]=wavelengthbd[i]+(double)centerblocker-lam0;
00887       blockerd[i]=blockerd[i]/100.0;
00888     }
00889   
00890   double frontwindowint[nlam];
00891   double blockerint[nlam];
00892   lininterp1f(frontwindowint,wavelengthd,frontwindowd,lam,ydefault2,nfront,nlam);   //Interpolation on the same wavelength grid
00893   lininterp1f(blockerint,wavelengthbd,blockerd,lam,ydefault2,nblocker,nlam);
00894   for(j=0;j<nlam;++j) blockerint[j]=blockerint[j]*frontwindowint[j];
00895 
00896   double lineprofile[nlam];
00897   double lineprofile2[referencenlam];
00898   int shiftlam;
00899   printf("Computing the HMI transmission profile\n");
00900 
00901 
00902  //OPENING OF BINARY FILES CONTAINING PHASES AND CONTRASTS OF NON-TUNABLE ELEMENTS
00903   int nelemPHASENT=4*nx2*nx2;
00904   float phaseNT[nelemPHASENT];    //phases of the non-tunable elements
00905   float contrastNT[nelemPHASENT]; //contrasts of the non-tunable elements
00906   int nread;
00907 
00908   printf("READ PHASES OF NON-TUNABLE ELEMENTS\n");
00909   if(nx2 == 256) fp = fopen(filephasesnontunable,"rb");    //in float, and 256x256. idl format: phase[256,256,4], so C format: phase[element][row][column]
00910   if(nx2 == 128) fp = fopen(filephasesnontunable128,"rb"); //in float, and 128x128.
00911   if(fp == NULL)
00912     {
00913       printf("CANNOT OPEN FILE OF NON-TUNABLE ELEMENT PHASES\n");
00914       exit(EXIT_FAILURE);
00915     }
00916   nread=fread(phaseNT,sizeof(float),nelemPHASENT,fp);
00917   fclose(fp);
00918   for(i=0;i<nelemPHASENT;++i) phaseNT[i] = phaseNT[i]*M_PI/180.; //convert phases from degrees to radians
00919 
00920   printf("READ CONTRASTS OF NON-TUNABLE ELEMENTS\n");
00921   if(nx2 == 256) fp = fopen(filecontrastsnontunable,"rb");    //in float, and 256x256. idl format: contrast[256,256,4], so C format: contrast[element][row][column]
00922   if(nx2 == 128) fp = fopen(filecontrastsnontunable128,"rb"); //in float, and 128x128.
00923   if(fp == NULL)
00924     {
00925       printf("CANNOT OPEN FILE OF NON-TUNABLE ELEMENT CONTRASTS\n");
00926       exit(EXIT_FAILURE);
00927     }
00928   nread=fread(contrastNT,sizeof(float),nelemPHASENT,fp);
00929   fclose(fp);
00930 
00931  //OPENING OF BINARY FILE CONTAINING CONTRASTS OF TUNABLE ELEMENTS
00932 
00933   int nelemCONTRASTT=3*nx2*nx2;
00934   float contrastT[nelemCONTRASTT]; //contrasts of the tunable elements
00935   printf("READ CONTRASTS OF TUNABLE ELEMENTS\n");
00936   if(nx2 == 256) fp = fopen(filecontraststunable,"rb");    //in float, and 256x256. idl format: contrast[256,256,3], so C format: contrast[element][row][column]
00937   if(nx2 == 128) fp = fopen(filecontraststunable128,"rb"); //in float, and 128x128. idl format
00938   if(fp == NULL)
00939     {
00940       printf("CANNOT OPEN FILE OF NON-TUNABLE ELEMENT PHASES\n");
00941       exit(EXIT_FAILURE);
00942     }
00943   nread=fread(contrastT,sizeof(float),nelemCONTRASTT,fp);
00944   fclose(fp);
00945 
00946 
00947   double lyot[nlam];
00948   double cmich[3];
00949   double inten[N];
00950   double vel1;
00951   double vel2;
00952   double filters[N][nlam];
00953   double values[3];
00954   int location;
00955   double phaseE2,phaseE3,phaseE4,phaseE5,contrastE2,contrastE3,contrastE4,contrastE5;
00956   double phaseNB,phaseWB,phaseE1,contrastNB,contrastWB,contrastE1;
00957   for(i=0;i<3;++i) cmich[i] = 2.0*M_PI/FSR[i];
00958 
00959   int maxshiftlam=round(vtest[ntest-1]*dlamdv/dlam);
00960           
00961   /***********************************************************************************************************/
00962   /*COMPUTE THE LOOK-UP TABLES FOR EACH PIXEL                                                                */
00963   /***********************************************************************************************************/
00964 
00965   printf("Computing the look-up tables (can take a few hours or so on n02, with 8 threads)\n");
00966 #pragma omp parallel for default(none) private(location,iii,i,j,k,values,lineprofile,shiftlam,inten,filters,f1c,f1s,vel1,f2c,f2s,vel2,row,column,FWHM,minimum,wavelength2,lineprofile2,lyot,phaseE2,phaseE3,phaseE4,phaseE5,contrastE2,contrastE3,contrastE4,contrastE5,contrastNB,contrastWB,contrastE1,phaseNB,phaseWB,phaseE1) shared(FSR,blockerint,cosi,sini,cos2i,sin2i,distance,nx2,ny2,ntest,nlam,lam,vtest,dlam,dlamdv,wavelength,N,cmich,phases,pv1,pv2,vel,axisout,nelement,ydefault,minimumCoeffs,FWHMCoeffs,WRONGDISTANCE,BUFFERDISTANCE,HCME1phase,HCMWBphase,HCMNBphase,phaseNT,contrastNT,contrastT,maxshiftlam,wavelengthref,referencenlam,shiftw,solarradiusmax,NOMINALSCALE,templineref)
00967   for(iii=0;iii<nx2*ny2;++iii)
00968     {
00969       row   =iii / nx2; //nx2= number of columns
00970       column=iii % nx2;
00971 
00972       //if(distance[iii] <= solarradiusmax*NOMINALSCALE)
00973       if(distance[iii] != WRONGDISTANCE)
00974     {
00975 
00976       if(column == (nx2/2)) printf("row %d column %d distance %f\n",row,column,distance[iii]);
00977 
00978 
00979       location  =column+row*nx2;
00980 
00981       //INTERPOLATION OF THE FeI LINEWIDTH AND LINEDEPTH AT DIFFERENT ANGULAR DISTANCES FROM DISK CENTER; BASED ON R. ULRICH'S LINES      
00982       FWHM=FWHMCoeffs[0]+FWHMCoeffs[1]*distance[iii];
00983       minimum=minimumCoeffs[0]+minimumCoeffs[1]*distance[iii];
00984       for(i=0;i<referencenlam;++i)
00985         {
00986           lineprofile2[i] = (1.0-(double)templineref[i])*minimum/(minimumCoeffs[0]+minimumCoeffs[1]); //scaling by the ratio of linedepths
00987           lineprofile2[i] =  1.0-lineprofile2[i];
00988           wavelength2[i]  = (double)wavelengthref[i]*FWHM/(FWHMCoeffs[0]+FWHMCoeffs[1]);
00989         }
00990       
00991       //INTERPOLATION OF THE FeI LINEWIDTH AND LINEDEPTH AT DIFFERENT ANGULAR DISTANCES FROM DISK CENTER; BASED ON GAUSSIAN PROFILES
00992 
00993       //FWHM=100.67102+0.015037016*distance[iii]-0.00010128197*distance[iii]*distance[iii]+3.1548385E-7*distance[iii]*distance[iii]*distance[iii]-3.7298102E-10*distance[iii]*distance[iii]*distance[iii]*distance[iii]+1.7275788E-13*distance[iii]*distance[iii]*distance[iii]*distance[iii]*distance[iii];// FWHM in mA
00994       //FWHM=FWHM/2.0/sqrt(log(2.0))/1000.; //we convert from FWHM in mA to sigma in A
00995       //minimum=0.515347+2.16653E-6*distance[iii]-2.76512E-7*distance[iii]*distance[iii]+2.92056E-10*distance[iii]*distance[iii]*distance[iii]-1.18263E-13*distance[iii]*distance[iii]*distance[iii]*distance[iii]; //linedepth for a continuum of 1
00996 
00997       /*distance[iii]=;
00998       FWHM=118.429+0.015037016*distance[iii]-0.00010128197*distance[iii]*distance[iii]+3.1548385E-7*distance[iii]*distance[iii]*distance[iii]-3.7298102E-10*distance[iii]*distance[iii]*distance[iii]*distance[iii]+1.7275788E-13*distance[iii]*distance[iii]*distance[iii]*distance[iii]*distance[iii];// FWHM in mA
00999       FWHM=FWHM/2.0/sqrt(log(2.0))/1000.; //we convert from FWHM in mA to sigma in A
01000       minimum=0.642710+2.16653E-6*distance[iii]-2.76512E-7*distance[iii]*distance[iii]+2.92056E-10*distance[iii]*distance[iii]*distance[iii]-1.18263E-13*distance[iii]*distance[iii]*distance[iii]*distance[iii]; //linedepth for a continuum of 1
01001 
01002 
01003 
01004       for(i=0;i<referencenlam;++i)
01005         {
01006           lineprofile2[i] = (1.0-(double)templineref[i])*minimum/minimumCoeffs; //scaling by the ratio of linedepths
01007           lineprofile2[i] =  1.0-lineprofile2[i];
01008           wavelength2[i]  = (double)wavelengthref[i]*FWHM/118.429;
01009           }*/
01010 
01011 
01012       //linear interpolation function (lininterpf1): for whatever reason, OpenMp does not work if I call linterpf1, therefore I have to copy this function here!!!
01013       for (i=0; i<nlam; i++)
01014         {
01015           if((lam[i] < (double)wavelength2[0]) || (lam[i] > (double)wavelength2[referencenlam-1])) lineprofile[i] = ydefault;
01016           else
01017         {   
01018           for(j=1; j<referencenlam; j++)
01019             {      
01020               if(lam[i]<=(double)wavelength2[j])
01021             {          
01022               lineprofile[i] = (lam[i]-(double)wavelength2[j-1]) / ((double)wavelength2[j]-(double)wavelength2[j-1]) * ((double)lineprofile2[j]-(double)lineprofile2[j-1]) + (double)lineprofile2[j-1];
01023               break;
01024             }
01025             }
01026         }
01027         }
01028 
01029 
01030       //for(k=0;k<nlam;++k) lineprofile[k]=1.0-phases[4+location*5]*minimum/minimumCoeffs*exp(-lam[k]*lam[k]/(phases[3+location*5]*FWHM/FWHMCoeffs)/(phases[3+location*5]*FWHM/FWHMCoeffs));
01031 
01032       //for(k=0;k<nlam;++k) lineprofile[k]=1.0-minimum*exp(-lam[k]*lam[k]/FWHM/FWHM);
01033 
01034 
01035       //NON-TUNABLE HMI TRANSMISSION PROFILE
01036       contrastNB=(double)contrastT[location];
01037       contrastWB=(double)contrastT[location+nx2*ny2];
01038       contrastE1=(double)contrastT[location+2*nx2*ny2];
01039       contrastE2=(double)contrastNT[location];
01040       contrastE3=(double)contrastNT[location+nx2*ny2];
01041       contrastE4=(double)contrastNT[location+2*nx2*ny2];
01042       contrastE5=(double)contrastNT[location+3*nx2*ny2];
01043       //phaseNB   =(double)phases[0+location*nelement];
01044       //phaseWB   =(double)phases[1+location*nelement];
01045       //phaseE1   =(double)phases[2+location*nelement];
01046       phaseNB   =(double)phases[0+location*5]*M_PI/180.;
01047       phaseWB   =(double)phases[1+location*5]*M_PI/180.;
01048       phaseE1   =(double)phases[2+location*5]*M_PI/180.;
01049       phaseE2   =(double)phaseNT[location];
01050       phaseE3   =(double)phaseNT[location+nx2*ny2];
01051       phaseE4   =(double)phaseNT[location+2*nx2*ny2];
01052       phaseE5   =(double)phaseNT[location+3*nx2*ny2];
01053 
01054 
01055       //iii = column, jjj = row
01056       //if(column == 100 && row == 50) printf("PHASES NT= %f %f %f %f %f %f %f %f %f %f %f\n",phaseE2,phaseE3,phaseE4,phaseE5,contrastNB,contrastWB,contrastE1,contrastE2,contrastE3,contrastE4,contrastE5);
01057       //exit(EXIT_FAILURE);
01058       
01059 
01060       //MODIFICATION: SHIFT OF NON-TUNABLE ELEMENTS PROFILES IN WAVELENGTH: (CAREFUL: FSR=2 pi/FSR)
01061       //phaseE2+=0.01432;//shiftw*FSR[3];
01062       //phaseE3+=0.055889;//shiftw*FSR[4];
01063       //phaseE4+=0.225455;//shiftw*FSR[5];
01064       //phaseE5+=0.900233;//shiftw*FSR[6];
01065 
01066       //phaseE2+=0.2;
01067       //phaseE3+=0.2;
01068       phaseE4+=0.4;
01069       phaseE5-=1.1;//0.900233;
01070 
01071 
01072       //CORRECTION OBS-CAL
01073       phaseNB+=1.573*M_PI/180.;
01074       phaseWB+=0.59*M_PI/180.;
01075       phaseE1+=-3.27*M_PI/180.;
01076       //phaseE2+=0.6*M_PI/180.;
01077       //phaseE3+=0.6*M_PI/180.;
01078       //phaseE4+=-1.1*M_PI/180.;
01079       //phaseE5+=-1.0*M_PI/180.;
01080 
01081 
01082       //MODIFICATION: SHIFT OF NON-TUNABLE ELEMENTS PROFILES IN WAVELENGTH:
01083       //phaseE2+=shiftw*FSR[3];
01084       //phaseE3+=shiftw*(FSR[4]/FSR[3])*FSR[4];
01085       //phaseE4+=shiftw*(FSR[5]/FSR[3])*FSR[5];
01086       //phaseE5+=shiftw*(FSR[6]/FSR[3])*FSR[6];
01087       
01088       for(j=0;j<nlam;++j) lyot[j]=blockerint[j]*(1.+contrastE2*cos(2.0*M_PI/FSR[3]*lam[j]+phaseE2))/2.*(1.+contrastE3*cos(2.0*M_PI/FSR[4]*lam[j]+phaseE3))/2.*(1.+contrastE4*cos(2.0*M_PI/FSR[5]*lam[j]+phaseE4))/2.*(1.+contrastE5*cos(2.0*M_PI/FSR[6]*lam[j]+phaseE5))/2.;
01089       
01090 
01091       //TUNABLE TRANSMISSION PROFILE
01092       for(i=0;i<N;i++) 
01093         {
01094           for(j=0;j<nlam;j++) filters[i][j] = lyot[j]*(1.+contrastNB*cos(cmich[0]*lam[j]+HCMNBphase[i]+phaseNB))/2.*(1.+contrastWB*cos(cmich[1]*lam[j]+HCMWBphase[i]+phaseWB))/2.*(1.+contrastE1*cos(cmich[2]*lam[j]-HCME1phase[i]+phaseE1))/2.;
01095         }
01096 
01097       //Computation of the norm of the filters
01098       //for(i=0;i<=N-1;++i) 
01099         //{
01100         //for(j=0;j<=nlam-1;++j) templine[j] = filters[i][j];
01101         //norm[i]=integration(lam,templine,nlam);
01102         //printf("Norm filter %d= %f\n",i,norm[i]);
01103         //}
01104       
01105       //look-up tables
01106 
01107 
01108 
01109       //START LOOP OVER THE INPUT VELOCITIES
01110       for(i=0;i<ntest;i++)
01111         {
01112 
01113           //BECAUSE dlam=dvtest*dlamdv, WE DO NOT NEED TO INTERPOLATE BUT WE CAN JUST SHIFT THE SOLAR LINE FOR DIFFERENT DOPPLER VELOCITIES (SLOW, NEED TO OPTIMIZE)
01114           shiftlam=round(vtest[i]*dlamdv/dlam);//SIGN CONVENTION: POSITIVE VELOCITIES CORRESPOND TO REDSHIFT (MOVEMENTS AWAY FROM OBSERVER)    
01115  
01116           for(j=0;j<N;j++)
01117         { 
01118           inten[j] = 0.;
01119           for(k=maxshiftlam;k<nlam-maxshiftlam;k++) inten[j] = inten[j]+filters[j][k]*lineprofile[k-shiftlam]; //shiftlam should correspond at the most to 154.44 mA (+7.5 km/s)
01120           //inten[j] = inten[j]/norm[j];
01121         }
01122 
01123 
01124           //First and Second Fourier coefficients
01125           f1c=0.0;
01126           f1s=0.0;
01127           f2c=0.0;
01128           f2s=0.0;
01129           for(j=0;j<N;j++) 
01130         {
01131           f1c += cosi[j] *inten[j];
01132           f1s += sini[j] *inten[j];
01133           f2c += cos2i[j]*inten[j];
01134           f2s += sin2i[j]*inten[j];
01135         }
01136 
01137           vel1 = atan2(-f1s,-f1c)*pv1/2.0/M_PI;
01138           vel[i+column*axisout[0]+row*axisout[0]*axisout[1]]       = (float)vel1;
01139           //vel[i+column*axisout[0]+row*axisout[0]*axisout[1]]       = (fmod((vel1-vtest[i]+10.5*pv1),pv1)-pv1/2.0+vtest[i]); //vel1 contains the "look-up" table
01140           vel2 = atan2(-f2s,-f2c)*pv2/2.0/M_PI;
01141           vel[i+ntest+column*axisout[0]+row*axisout[0]*axisout[1]] = (float)(fmod((vel2-vel1+10.5*pv2),pv2)-pv2/2.0+vel1);           //vel2 contains the "look-up" table
01142 
01143         }//for(i=0;i<ntest;i++)
01144       
01145     } //if(distance[iii] <= 1.0 && distance[iii] >= 0.0)   
01146     }
01147 
01148 
01149   //SUPPORT FOR LARGER CROP RADII
01150   double mean[ntest],mean2[ntest];
01151   for(row=0;row<nx2;++row)
01152     {
01153       if(row < nx2/2)
01154     {
01155       for(i=0;i<ntest;i++)
01156         {
01157           mean[i] =vel[i+nx2/2*axisout[0]+1*axisout[0]*axisout[1]];
01158           mean2[i]=vel[i+ntest+nx2/2*axisout[0]+1*axisout[0]*axisout[1]];
01159         }
01160     }
01161       else
01162     {
01163       for(i=0;i<ntest;i++)
01164         {
01165           mean[i] =vel[i+nx2/2*axisout[0]+(nx2-2)*axisout[0]*axisout[1]];
01166           mean2[i]=vel[i+ntest+nx2/2*axisout[0]+(nx2-2)*axisout[0]*axisout[1]];
01167         }
01168     }
01169 
01170       for(column=nx2/2;column>=0;--column)
01171     {
01172       location  =column+row*nx2;
01173       if(distance[location] != WRONGDISTANCE)
01174         {
01175           for(i=0;i<ntest;i++)
01176         {
01177           mean[i] =vel[i+column*axisout[0]+row*axisout[0]*axisout[1]];
01178           mean2[i]=vel[i+ntest+column*axisout[0]+row*axisout[0]*axisout[1]];
01179         }
01180         }
01181       if(distance[location] == WRONGDISTANCE)
01182         {
01183           for(i=0;i<ntest;i++)
01184         {
01185           vel[i+column*axisout[0]+row*axisout[0]*axisout[1]]       = mean[i];
01186           vel[i+ntest+column*axisout[0]+row*axisout[0]*axisout[1]] = mean2[i];
01187         }
01188         }
01189     }
01190 
01191       for(column=nx2/2;column<=nx2-1;++column)
01192     {
01193       location  =column+row*nx2;
01194       if(distance[location] != WRONGDISTANCE)
01195         {
01196           for(i=0;i<ntest;i++)
01197         {
01198           mean[i] =vel[i+column*axisout[0]+row*axisout[0]*axisout[1]];
01199           mean2[i]=vel[i+ntest+column*axisout[0]+row*axisout[0]*axisout[1]];
01200         }
01201         }
01202       if(distance[location] == WRONGDISTANCE)
01203         {
01204           for(i=0;i<ntest;i++)
01205         {
01206           vel[i+column*axisout[0]+row*axisout[0]*axisout[1]]       = mean[i];
01207           vel[i+ntest+column*axisout[0]+row*axisout[0]*axisout[1]] = mean2[i];
01208         }
01209         }
01210     }
01211     }
01212 
01213   /***********************************************************************************************************/
01214   /* SAVE THE RESULTS                                                                                        */
01215   /***********************************************************************************************************/
01216 
01217   //create a record in the series dsout (Look-up tables)
01218   DRMS_RecordSet_t *dataout = NULL;
01219   DRMS_Record_t  *recout    = NULL;
01220   DRMS_Segment_t *segout    = NULL;
01221 
01222   dataout = drms_create_records(drms_env,1,dsout,DRMS_PERMANENT,&status);
01223   
01224   if (status != DRMS_SUCCESS)
01225     {
01226       printf("Could not create a record for the lookup tables\n");
01227       exit(EXIT_FAILURE);
01228     }
01229   if (status == DRMS_SUCCESS)
01230     {     
01231       printf("Writing a record on the DRMS for the lookup tables\n");
01232       recout = dataout->records[0];
01233  
01234       //WRITE KEYWORDS
01235       status = drms_setkey_int(recout,keyname,keyvalue);
01236       status = drms_setkey_time(recout,keyname2,interntime);
01237       status = drms_setkey_int(recout,keyname3,keyvalue3);
01238       status = drms_setkey_time(recout,keyname4,interntime2);
01239       status = drms_setkey_int(recout,keyname6,keyvalue5);
01240       status = drms_setkey_time(recout,keyname5,interntime3);
01241       status = drms_setkey_int(recout,keyname7,camera);//HCAMID
01242       status = drms_setkey_double(recout,CROTA2,crota2);  // added for ROLL 2018
01243       status = drms_setkey_double(recout,OBS_VR,obs_vr);  // added for ROLL 2018
01244       char  DATEOBS[256];
01245       sprint_time(DATEOBS,CURRENT_SYSTEM_TIME,"UTC",1);
01246       status = drms_setkey_string(recout,DATEs,DATEOBS);           //DATE AT WHICH THE FILE WAS CREATED
01247       status = drms_setkey_string(recout,COMMENTS,COMMENT);        //COMMENT ON THE DATA 
01248       status = drms_setkey_float(recout,FSRNBs,FSR[0]);
01249       status = drms_setkey_float(recout,FSRWBs,FSR[1]);
01250       status = drms_setkey_float(recout,FSRE1s,FSR[2]);
01251       status = drms_setkey_float(recout,FSRE2s,FSR[3]);
01252       status = drms_setkey_float(recout,FSRE3s,FSR[4]);
01253       status = drms_setkey_float(recout,FSRE4s,FSR[5]);
01254       status = drms_setkey_float(recout,FSRE5s,FSR[6]);
01255       status = drms_setkey_float(recout,CBLOCKERs,centerblocker);
01256 
01257       //FOLLOWING LINES ARE FOR THE ACCIDENTAL REBOOT OF APRIL 2013. IT SHOULD BE COMMENTED FOR ANY LOOK-UP TABLE THAT IS NOT AT THAT TIME
01258       //FOR FSN_REC=55032512 CORRESPONDING TO HMI ACCIDENTAL REBOOT OF END OF APRIL 2013, I FORCE N=7 WHEN CREATING THE LOOKUP TABLE
01259       /*HCMNB=81;
01260       HCMWB=52;
01261       HCME1=37;
01262       HCMPOL=0;
01263       N=7;*/
01264 
01265       status = drms_setkey_int(recout,HCMNBs,HCMNB);
01266       status = drms_setkey_int(recout,HCMWBs,HCMWB);
01267       status = drms_setkey_int(recout,HCME1s,HCME1);
01268       status = drms_setkey_int(recout,HCMPOLs,HCMPOL);
01269       status = drms_setkey_int(recout,Ns,N);
01270       status = drms_setkey_float(recout,RADIUS,solarradiusmax);
01271 
01272       //WRITE DATA SEGMENT
01273       segout = drms_segment_lookupnum(recout, 0);
01274       drms_segment_write(segout, arrout, 0);           //write the file containing look-up tables
01275 
01276       //CLOSE RECORDS
01277       drms_close_records(dataout, DRMS_INSERT_RECORD); //insert the record in DRMS
01278     }
01279   
01280   drms_free_array(arrin);
01281   drms_free_array(arrout);
01282   free(cosi);
01283   free(sini);
01284   free(cos2i);
01285   free(sin2i);
01286   free(angle);
01287   free(HCME1phase);
01288   free(HCMWBphase);
01289   free(HCMNBphase);
01290   drms_close_records(data, DRMS_FREE_RECORD); //insert the record in DRMS
01291   
01292   return(0);
01293   
01294 }
proj/lev1.5_hmi/apps/HMI_lookup.c
