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

00001 // USING A VOIGT PROFILE FOR THE Fe I LINE
00002 
00003 /*-------------------------------------------------------------------------------------------------------*/
00004 
00005 /*  Modified from JSOC/proj/lev1.5_hmi/apps/phasemaps_voight.c
00006 
00007 /*-------------------------------------------------------------------------------------------------------*/
00008 
00009 /*-------------------------------------------------------------------------------------------------------
00010  * Changes needed for OBS_MODE  -- Phil Jan 2016
00011  * OBS_MODE detunes use the same sequence, 27 exposures with a dark at each end.
00012  * In OBS_MODE the WCS keywords are populated.
00013  * Due to CRPIX2 jitter of 0.75 pixel and about half that for CRPIX1 the image must
00014  * be smoothed perhaps as well as binning to eliminate high gradients in LOS velocity
00015  * Changes:
00016  *   0.  Use HCFTID to determine CAL_MODE or OBS_MODE
00017  *   1.  Must compute LOS velocity for each location
00018  *   2.  Must set solarradiusmax to R_SUN in right units...
00019  *-------------------------------------------------------------------------------------------------------
00020  */
00021 /* XXXXXX Changes made in Jan 2016 by Phil
00022  * 1. changed comment for HCAMID
00023  * 2. Changes all instances in DoIt of 'exit(' to 'return('
00024  * 3. Add code from sunrotation.c for orbit removal and average sun doppler removal
00025  */
00026 
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include <math.h>
00030 #include <gsl/gsl_sort_double.h>//GSL = Gnu Scientifc Library
00031 #include <gsl/gsl_blas.h>       //GNU implementation of the BLAS routines
00032 #include <gsl/gsl_linalg.h>     //SVD algorithm from the GNU Scientific Library
00033 #include <gsl/gsl_statistics.h> //mean and sigma
00034 
00035 #include "mdi.h"           //contains definitions for some HMI filter functions
00036 #include "mdi_params.h"      //contains definitions for some HMI filter parameters
00037 #include "mdi.c"           //contains some HMI filter functions
00038 
00039 double Ni_Line   =  LAMBDA_Ni; //Ni I line central wavelength in vacuum at rest, 6767.78 A 
00040 
00041 #include <omp.h>                //Open MP header
00042 #include <fresize.h>            //from Jesper: to rebin the 1024x1024 images
00043 #include "interpol_code.h"      //from Richard, for de-rotation and gap-filling
00044 
00045 char *module_name    = "MDI_phasemaps";   //name of the module
00046 
00047 #define kRecSetIn      "input_series" //names of the arguments of the module
00048 #define kDSOut         "phasemap_series"
00049 #define kreduced       "reduced"      //maps in 64x64 or 256x256?             
00050 #define minval(x,y) (((x) < (y)) ? (x) : (y))
00051 #define maxval(x,y) (((x) < (y)) ? (y) : (x))
00052 
00053 #define Deg2Rad  (M_PI/180.0)
00054 
00055 //arguments of the module
00056 ModuleArgs_t module_args[] =        
00057 {
00058      {ARG_STRING, "in", "" ,  "Input data series."},
00059      {ARG_STRING, "out",    "" ,  "Phase Maps series."},
00060      {ARG_INT   , kreduced,  "0",  "64x64 (1), 32x32 (2), 128x128 (3), or standard (256x256) resolution (0)?"},
00061      {ARG_DOUBLE , "FSRNB", "-9999.0" , "FSR NB"},
00062      {ARG_DOUBLE , "FSRWB", "-9999.0" , "FSR WB"},
00063      {ARG_FLOAT  , "center", "" , "center of blocker filter"},
00064      {ARG_FLOAT  , "shift", "" , "wavelength shift of the non-tunable profile"},
00065      {ARG_DOUBLE , "thresh", "", "threshold"},
00066      {ARG_END}
00067 };
00068 
00069 // function to get LOS velocity of observed point
00070 double Solar_velocity(double x,double y,double rsun2,double OBS_B0,double CROTA2, double CDELT, double obs_vr, double obs_vw, double obs_vn);
00071 int img2sphere (double x, double y, double ang_r, double latc, double lonc,
00072     double pa, double *p_rho, double *p_lat, double *p_lon, double *p_sinlat,
00073     double *p_coslat, double *p_sig, double *p_mu, double *p_chi);
00074 // rebin function that deals with missing data.
00075 int rebinArraySF(float *outData, float *inData, int outNx, int inNx, float rsun, float x0, float y0, float fillvalue);
00076 
00077 /*-----------------------------------------------------------------------------------------------------*/
00078 /* Function to perform linear interpolation                                                            */
00079 /* found on the internet                                                                               */
00080 /* returns the values yinterp at points x of 1D function yv (at the points xv)                         */
00081 /*-----------------------------------------------------------------------------------------------------*/
00082 
00083 void lininterp1f(double *yinterp, double *xv, double *yv, double *x, double ydefault, int m, int minterp)
00084 {
00085     int i, j; 
00086     int nrowsinterp, nrowsdata;
00087     nrowsinterp = minterp;
00088     nrowsdata = m;
00089     for (i=0; i<nrowsinterp; i++)
00090       {
00091     if((x[i] < xv[0]) || (x[i] > xv[nrowsdata-1])) yinterp[i] = ydefault;
00092     else
00093       {   
00094         for(j=1; j<nrowsdata; j++)
00095           {      
00096         if(x[i]<=xv[j])
00097           {        
00098             yinterp[i] = (x[i]-xv[j-1]) / (xv[j]-xv[j-1]) * (yv[j]-yv[j-1]) + yv[j-1];
00099             break;
00100           }
00101           }
00102       }
00103       }
00104 } 
00105 
00106 
00107 /*-----------------------------------------------------------------------------------------*/
00108 /*                                                     */
00109 /* FUNCTION TO PRODUCE A LOOK-UP TABLE CORRESPONDING TO THE PHASES                     */
00110 /*                                                     */
00111 /* the convention for the filter names is (with 4 filters for Doppler):                    */
00112 /* F1 is at -113.4 mA   M1 @ -27= -0.6 FSR1, M2 @ -13=-0.2889 FSR2                         */
00113 /* F2 is at  -37.8 mA   M1 @ -9 = -0.2 FSR1, M2 @ -4 =-0.0889 FSR2                         */
00114 /* F3 is at   37.8 mA   M1 @  9 =  0.2 FSR1, M2 @  5 = 0.1111 FSR2                         */
00115 /* F4 is at  113.4 mA   M1 @  27=  0.6 FSR1, M2 @ 14 = 0.3111 FSR2                         */
00116 /* Fc is at -184.3 mA   M1 @ -45= -1.0 FSR1, M2 @-22 =-0.4889 FSR2                         */
00117 /*                                                                                         */
00118 /* the relative phases of NB and WB increase when the steps of the HCMs                    */
00119 /* increase                                                                                */
00120 /*-----------------------------------------------------------------------------------------*/
00121 
00122 void cotunetable(double NBphase,double WBphase,int table[2][20])
00123 {
00124 
00125   int HCMNB, HCMWB, i;
00126 
00127   HCMNB    = (int)round(-NBphase/6.0)+60;            //for convention Intensity = 0.5*[1+cos(lam+phase)]
00128   HCMWB    = (int)round(-WBphase/6.0)+60;
00129 
00130 
00131   table[0][0]=-30; table[1][0]=0;
00132   table[0][1]=-27; table[1][1]=-6;
00133   table[0][2]=-24; table[1][2]=-12;
00134   table[0][3]=-21; table[1][3]=-18;
00135   table[0][4]=-18; table[1][4]=-24;
00136   table[0][5]=-15; table[1][5]=-30;    //I5 for N=6 tuning positions
00137   table[0][6]=-12; table[1][6]=24;   //I4 for N=5
00138   table[0][7]=-9;  table[1][7]=18;   //I4 for N=6
00139   table[0][8]=-6;  table[1][8]=12;   //I3 for N=5
00140   table[0][9]=-3;  table[1][9]=6;    //I3 for N=6
00141   table[0][10]=0;  table[1][10]=0;   //I2 for N=5
00142   table[0][11]=3;  table[1][11]=-6; //I2 for N=6
00143   table[0][12]=6;  table[1][12]=-12;//I1 for N=5
00144   table[0][13]=9;  table[1][13]=-18; //I1 for N=6
00145   table[0][14]=12; table[1][14]=-24;//I0 for N=5
00146   table[0][15]=15; table[1][15]=-30; //I0 for N=6
00147   table[0][16]=18; table[1][16]=24; 
00148   table[0][17]=21; table[1][17]=18;
00149   table[0][18]=24; table[1][18]=12;
00150   table[0][19]=27; table[1][19]=6;
00151 
00152 
00153   for(i=0;i<20;++i)
00154     {
00155       table[0][i] = table[0][i] + HCME1;
00156       table[1][i] = table[1][i] + HCMWB;
00157     }
00158 
00159 }
00160   
00161 
00162 /*------------------------------------------------------------------------------------------------------*/
00163 /*                                                                                                      */
00164 /*  MAIN PROGRAM                                                                                        */
00165 /*                                                                                                      */
00166 /*------------------------------------------------------------------------------------------------------*/
00167 
00168 
00169 #define UNK_MODE 0
00170 #define CAL_MODE 1
00171 #define OBS_MODE 2
00172 
00173 #define DIE(msg) {fprintf(stderr,msg); return(EXIT_FAILURE); }
00174 #define DIE1(msg,val) {fprintf(stderr,msg,val); return(EXIT_FAILURE); }
00175 
00176 int DoIt(void)
00177   {
00178   const char *inRecQuery     = cmdparams_get_str(&cmdparams, "in", NULL); //cmdparams is defined in jsoc_main.h
00179   const char *dsout          = cmdparams_get_str(&cmdparams, "out" ,NULL);    //series name of the output phase maps series
00180   int reduced          = cmdparams_get_int(&cmdparams  , kreduced ,NULL);  //if reduced=1 then we only compute 64x64 maps (nx2=64), =2 it's 32x32, =3 128x128, instead of the usual 256x256 (reduced =0)
00181 
00182   // Parameters defined in HMIparam.h with values for MDI
00183   float centerblocker2 =  cmdparams_get_float(&cmdparams,"center", NULL);
00184   float shiftw         = cmdparams_get_float(&cmdparams, "shift" , NULL);
00185   double FSR[2] = {FSR1,FSR2};
00186   int nx2 = 128;        // phasemap dimensions
00187   int ny2 = 128;        // phasemap dimensions
00188   int nelement = 2;     // number tunable elements
00189   int ntest = 801;      // number of test velocities
00190   double dvtest = 10.0      // m/s, so -4000:4000 m/s 
00191   double vtestmax = 4000.0
00192   int nfront = 201;     // number of wavelengths in front window transmission profile
00193   int nblocker = 401;       // number of wavelengths in blocker and Lyot profile, +- 1A so 5 mA resolution.
00194   double centerblocker = LCEN   // Lyot center wavelength WRT Ni line 0 velocity.
00195   double dlamdv = Ni_Line/C;    // 2.2574e-5 Angstroms per m/s
00196   int numberoffiltergrams_detune = 20 ; // 9 for detune, 1 FC, 10 for Dop sequence.
00197   double solar_radius = 696000000.0;
00198   double Astrounit = 149597870691.0;
00199   float phaseguess[3] = {272.0, 24.0};  // initial guess for cotune phase, assuming motor positions 79, 48
00200   double depth0 = 0.64;     // Ni line depth (I0 = 0.4)
00201   double widtho = 0.116;    // Ni Line FWHM Angstroms, (Norton et al)
00202   double thresh = 60000;    // guess for continuum
00203   float MINRSUN = 0.0;
00204   float MAXRSUN = 512;
00205   float XYMIN = 0;
00206   float XYMAX = 1024;
00207   double NOMINALSCALE = 2.0;    // scale for calmode, arcsec/pixel
00208 
00209   // Defaults for the following 8 params are in HMIparam.h
00210   double try;
00211   try = cmdparams_get_double(&cmdparams,"FSRNB" , NULL);
00212     if (try != -9999.0) FSR[0] = try;
00213   try = cmdparams_get_double(&cmdparams,"FSRWB" , NULL);
00214     if (try != -9999.0) FSR[1] = try;
00215   try = cmdparams_get_double(&cmdparams,"thresh", NULL);
00216     if (try != -9999.0) thresh = try;
00217   printf("FSR[0]=%f, FSR[1] = %f\n",FSR[0],FSR[1]);
00218 
00219   char COMMENT[256];
00220   
00221   char *COMMENTS= "COMMENT";
00222 
00223   printf("FSR PARAMETERS = %f %f %f %f %f\n",FSR[0],centerblocker2,shiftw,thresh);
00224   printf("COMMENT: %s\n",COMMENT);
00225 
00226   if(reduced == 1)
00227     {
00228       nx2=64;
00229       ny2=64;
00230     }
00231   if(reduced == 2)
00232     {
00233       nx2=32;
00234       ny2=32;
00235     }
00236   if(reduced == 3)
00237     {
00238       nx2=128;
00239       ny2=128;
00240     }
00241 
00242   int  status       = DRMS_SUCCESS; 
00243   int  error        = 0;
00244   int  nx           = 1024;                                           //number of columns of input lev 1 filtergrams
00245   int  ny           = 1024;                                           //number of rows of input lev 1 filtergrams
00246   int  Nelem        = nx*ny;
00247   int  ratio;                                                         //ratio nx/nx2 (MUST BE AN EVEN NUMBER)
00248   double iterstop   = 1.e-7;                                //criterion to stop the iterations: we stop when the maximum relative change in the parameters is less than iterstop
00249 
00250 #define nseq  9                                                      //number of wavelength positions in the detune sequences (EXLUDING THE DARK FRAMES)
00251 #define maxsteps 155                                                //maximum steps allowed for the iterative Least-Squares algorithm
00252 #define nparam 6                                                      //number of parameters we fit for
00253 
00254   int nimages       = nseq+1;                                         //number of images in the detune sequence
00255   int indeximage[nimages];
00256   double factor     = 1000.;                                          //to make sure all the fitted parameters are, roughly, of the same order of magnitude
00257   double dpi        = 2.0*M_PI;
00258   int i,j,k,iii,jjj;                                                  //variables for the loops
00259   float distance2[nx*ny];
00260   float distance[ny2*nx2];                                            //format: [row][column] (in idl notation: distance[nx2,ny2])
00261                                                                       //nx2 is a number of column, ny2 is a number of rows
00262   double Inten[ny2][nx2][nseq];
00263   float solarradiusmax;                                               //maximum radius at which the phases will be computed
00264 
00265   //FOR DO_INTERPOLATE()
00266 #define TABLEPATH "/home/jsoc/cvs/Development/JSOC/proj/lev1.5_hmi/libs/lev15/"
00267 
00268   struct init_files initfiles;
00269   char MDISeriesTemp[256];
00270 
00271   //WAVELENGTH GRID PARAMETERS
00272 
00273 #define  nlam 20002                                                    //MUST BE AN EVEN NUMBER
00274 double dlam       = 0.0005;
00275 
00276   double lam[nlam];
00277   for(i=0;i<nlam;i++) lam[i] = ((double)i-((double)nlam-1.0)/2.0)*dlam;//symmetrical around 0 but NO 0.
00278  
00279   double blockerint[nlam];
00280   double templine[nlam];
00281   double line[nlam];
00282   double linea[nlam];
00283   double lineb[nlam];
00284   double linec[nlam];
00285   double lined[nlam];
00286   double gaussian[nlam];
00287   double l;
00288   double a=0.03; //damping for Voigt profile (value from Foy, 1972. Should be 0.0315 actually)
00289 
00290   double profilef[nlam];
00291   double profileg[nlam];
00292   double dlinedI0[nlam];
00293   double dlinedIc[nlam];
00294   double dlinedw[nlam];
00295   double err[maxsteps];
00296   double history[maxsteps][nparam];
00297   double ydefault    = 0.;
00298   int    converg;
00299   int    compteur;
00300   double dIntendI0[nseq];
00301   double dIntendw[nseq];
00302   double dIntendIc[nseq];
00303   double dIntendPhi0[nseq];
00304   double dIntendPhi1[nseq];
00305   double Icg[nx2][nx2];
00306   double fwidthg[nx2][nx2];
00307   double fdepthg[nx2][nx2];
00308   double residual;
00309   double jose;
00310   thresh=thresh/factor;   //we will divide the intensities by factor, thresh is defined in HMIparam.h
00311   size_t minimum=0;
00312   FILE *fp=NULL;
00313   int row, column;
00314   int indexref;
00315 
00316   int nthreads=omp_get_num_procs();         //number of threads supported by the machine where the code is running
00317   if (nthreads > 8) nthreads = 8;
00318   omp_set_num_threads(nthreads);            //set the number of threads to the maximum value
00319   printf("number of threads for OpenMP = %d\n",nthreads);
00320 
00321   //TEST THE REBINNING RATIO
00322   if( (nx % nx2) == 0)  ratio =nx/nx2;                                 //typically ratio=16, because nx2=256 and nx=1024
00323   else
00324     {
00325       printf("Error: nx = %d must be a multiple of nx2 = %d\n",nx,nx2);
00326       DIE("");
00327     }
00328 
00329   //for(i=0;i<nlam;++i) gaussian[i]=0.015*exp(-(lam[i]+0.225)*(lam[i]+0.225)/0.2/0.2)-0.004*exp(-(lam[i]-0.150)*(lam[i]-0.150)/0.22/0.22); //GAUSSIAN LINE ADDED TO VOIGT PROFILE TO SIMULATE ASYMMETRY
00330   //for(i=0;i<nlam;++i) gaussian[i]=0.015*exp(-(lam[i]+0.225)*(lam[i]+0.225)/0.2/0.2);
00331   //for(i=0;i<nlam;++i) gaussian[i]=0.0;
00332 
00333   /***********************************************************************************************************/
00334   /*DRMS ACCESS                                                                                              */
00335   /***********************************************************************************************************/
00336 
00337   printf("QUERY = %s\n",inRecQuery);
00338   DRMS_RecordSet_t *data = drms_open_records(drms_env,inRecQuery,&status);//open the records from the input series
00339 
00340   if (status == DRMS_SUCCESS && data != NULL && data->n > 0)
00341     {
00342       int nRecs = data->n;                                               //number of records in the input series 
00343       printf("Number of filtergrams to be read= %d \n",nRecs);           //check if the number of records is appropriate
00344       if (nRecs != numberoffiltergrams_detune)  
00345       DIE1("Problem: program requires %d filtergrams to produce the phase maps\n",numberoffiltergrams_detune);
00346 
00347       //the format of the detune9 sequence is assumed to be front camera, then side camera, then front camera, and so on...
00348       for(i=0;i<nimages;++i) indeximage[i]=i
00349 
00350       DRMS_Record_t *rec[nimages];                                       //we only work with 1 camera (9 detune positions)
00351       DRMS_Record_t *rec2[2]     ;                                       //records for the first and last filtergrams of the detune sequence
00352       for(i=0;i<nimages;++i) rec[i] = data->records[indeximage[i]];      //records that will be used to compute the phase maps
00353       rec2[0] = data->records[0];                                        //records that will be used to set the T_REC and FSN_REC keywords for the phasemaps
00354       rec2[1] = data->records[nimages-1];
00355 
00356       DRMS_Segment_t *segin  = NULL;
00357       int  Frecnum          = 0;
00358       int  Lrecnum         = 0;
00359       int  Mrecnum         = 0;  // average recnum
00360       DRMS_Type_t type       = DRMS_TYPE_FLOAT;
00361       DRMS_Type_t typeEr     = DRMS_TYPE_CHAR;
00362 
00363       const char *HCMNB      = "MTM1";                            //keyword for the position of the HCM of NB Michelson
00364       const char *HCMWB      = "MTM2";                            //keyword for the position of the HCM of WB Michelson
00365       const char *CRPIX1     = "CRPIX1";                             //x-axis location of solar disk center (IN OBS_MODE) in pixels, starting at 1
00366       const char *CRPIX2     = "CRPIX2";                             //y-axis location of solar disk center in pixels
00367       const char *SCALE      = "CDELT1";                             //image scale in the x direction
00368       const char *VR         = "OBS_VR";                             //SDO orbital velocity away from Sun m/s 
00369       const char *VW         = "OBS_VW";                             //SDO orbital velocity in direction of increasing HGLN_OBS 
00370       const char *VN         = "OBS_VN";                             //SDO orbital velocity in direction of increasing HGLT_OBS
00371       const char *RADIUS     = "R_SUN";                              //Image radius in pixels
00372       const char *EXPTIME    = "EXPTIME";                            //mean shutter open time
00373       const char *HCMNB0     = "HCMNB";
00374       const char *HCMWB0     = "HCMWB";
00375       const char *NXs        = "NX";
00376       const char *FSRNBs     = "FSRNB";
00377       const char *FSRWBs     = "FSRWB";
00378       const char *DATEs      = "DATE";
00379 
00380       int Mode               = UNK_MODE;
00381       TIME   interntime;
00382       TIME   interntime2;
00383       TIME   interntime3;
00384       float  X0[nseq],Y0[nseq];                                      //location of Sun's center in pixels (starts 0.0)
00385       float  CDELT1[nseq];                                           //image scale in the x direction
00386       float  RSUN[nseq];
00387       int    FSN,NBADPERM,imcfg;
00388       double VELOCITY[nseq];                                         //Sun-SDO radial velocity
00389       double EXPOSURE[nseq];
00390       double OBS_VR, OBS_VW, OBS_VN, lam0, dlam_LOS[nx2*ny2];
00391 
00392       //initialization of Richard's code for de-rotation and gapfilling, and Jesper's code for rebinning
00393       //************************************************************************************************
00394       
00395       struct  initial const_param;                                   //structure containing the parameters for Richard's functions
00396  
00397       char dpath[]="/home/jsoc/cvs/Development/JSOC/";
00398       status = initialize_interpol(&const_param,&initfiles,nx,ny,dpath);
00399       if(status != 0)
00400       DIE("Error: could not initialize the gapfilling routine\n");
00401       
00402       Mask = (unsigned char *)malloc(Nelem*sizeof(unsigned char));
00403       if(Mask == NULL)
00404       DIE("Error: cannot allocate memory for Mask\n");
00405 
00406       struct fresize_struct fresizes;
00407       init_fresize_bin(&fresizes,ratio);
00408 
00409       
00410       /***********************************************************************************************************/
00411       /*READING, GAPFILLING, AND REBINNING THE FILTERGRAMS                                                       */
00412       /***********************************************************************************************************/
00413 
00414       int tuningint[nseq][2]; //HCM positions of the tuning waveplates
00415       int axisout22[2] = {nx,ny}; //column, row
00416  
00417       float  tempframe[ny2*nx2];
00418       double tuning[nseq][];
00419       float  *arrinL = NULL;       
00420       char  *ierror  = NULL;
00421 
00422       DRMS_Array_t *arrin     = NULL;
00423       DRMS_Array_t *Ierror    = NULL;
00424  
00425       Ierror = drms_array_create(typeEr,2,axisout22,NULL,&status);
00426       if(status != DRMS_SUCCESS || Ierror == NULL)
00427       DIE("Error: could not create Ierror\n");
00428 
00429 
00430       //LOOP OVER ALL THE FILTERGRAMS OF THE DETUNE SEQUENCE
00431       printf("READING, GAPFILLING, AND REBINNING THE FILTERGRAMS in %dx%d\n",nx2,ny2);
00432       for(i=0;i<nseq;i++)                                            //read and rebin the filtergrams
00433     {
00434       //READING IMAGE
00435       segin     = drms_segment_lookupnum(rec[i],0);            //assume first record is first detune filtergram
00436       arrin     = drms_segment_read(segin,type,&status);
00437       if(status != DRMS_SUCCESS)
00438           DIE("Error: unable to read a data segment\n");
00439       arrinL = arrin->data;
00440 
00441       //READING SOME KEYWORDS
00442       FSN = drms_getkey_int(rec[i],"recnum",&status);
00443       printf("FSN IMAGE = %d\n",FSN);
00444       NBADPERM =  drms_getkey_int(rec[i],"NBADPERM",&status);
00445       if(status != DRMS_SUCCESS) NBADPERM=-1;
00446 
00447       tuningint[i][0] = drms_getkey_int(rec[i],"RLC_MTM1" ,&status); //tuning positions of the HCM of NB
00448       if(status != DRMS_SUCCESS)
00449           DIE1("Error: unable to read the keyword %s\n","MDI_MTM1");
00450       else
00451               printf(" %d ",tuningint[i][0]);
00452       tuningint[i][1] = drms_getkey_int(rec[i],"RLC_MTM2",&status); //tuning positions of the HCM of WB
00453       if(status != DRMS_SUCCESS)
00454           DIE1("Error: unable to read the keyword %s\n","MDI_MTM2");
00455       else
00456               printf(" %d ",tuningint[i][1]);
00457 
00458           DPC = drms_getkey_int(rec[i], "DPC", &status);
00459           if (i==0)
00460             {
00461         Mode = (DPC == 0x404e4400 ? CAL_MODE : OBS_MODE);
00462             printf(" %0X ",DPC);
00463             }
00464           if (Mode == OBS_MODE)
00465             {
00466         X0[i]         = drms_getkey_float(rec[i],"MDI_X0",&status);
00467         if(status != DRMS_SUCCESS || isnan(X0[i]) || X0[i] == -1)
00468            fprintf(stderr,"Error: unable to read the keyword %s\n","MDI_X0");
00469             else
00470            printf(" %f ",X0[i]);
00471         Y0[i]         = drms_getkey_float(rec[i],CRPIX2,&status);
00472         if(status != DRMS_SUCCESS || isnan(Y0[i]) || Y0[i] == -1)
00473             fprintf(stderr,"Error: unable to read the keyword %s\n",CRPIX2);
00474             else
00475            printf(" %f ",Y0[i]);
00476             }
00477           else
00478             {
00479         X0[i]= 511.5;
00480         Y0[i]= 511.5;
00481             }
00482 
00483       CDELT1[i]  = drms_getkey_float(rec[i],SCALE,&status); //I don't think I really need CDELT...
00484       if(status != DRMS_SUCCESS || isnan(CDELT1[i]))
00485         {
00486         fprintf(stderr,"Warning CDELT1 not found, use 0.5\n");
00487         CDELT1[i]=0.5;
00488         }
00489       printf(" %f ",CDELT1[i]);
00490 
00491           if (Mode == OBS_MODE)
00492             {
00493         RSUN[i]  = drms_getkey_float(rec[i],RADIUS,&status);
00494         if(status != DRMS_SUCCESS || isnan(RSUN[i]))
00495               DIE1("Error: unable to read the keyword %s\n",RADIUS);
00496             }
00497           else
00498          RSUN[i]=511.5;
00499       printf(" %f ",RSUN[i]);
00500 
00501       VELOCITY[i]= drms_getkey_double(rec[i],VR,&status);
00502       if(status != DRMS_SUCCESS || isnan(VELOCITY[i]))
00503           DIE1("Error: unable to read the keyword %s\n",VR);
00504           printf(" %f ",VELOCITY[i]);
00505 
00506       EXPOSURE[i]= drms_getkey_double(rec[i],EXPTIME,&status);
00507       if(status != DRMS_SUCCESS || isnan(EXPOSURE[i]))
00508           DIE1("Error: unable to read the keyword %s\n",EXPTIME);
00509        printf(" %f ",EXPOSURE[i]);
00510 
00511         // Rebin the input array
00512         if (Mode == OBS_MODE)
00513           {
00514           // Need to rebin allowing missing data
00515           rebinArraySF(tempframe, arrinL, nx2, nx, RSUN[i], X0[i], Y0[i], 0.0);
00516           }
00517         else
00518           {
00519           //REBINNING, USING THE FUNCTION OF JESPER
00520           fresize(&fresizes,arrinL,tempframe,nx,ny,nx,nx2,ny2,nx2,0,0,0.0f);
00521           }
00522 
00523         //put the image contained in tempframe into the Inten array [IN DOUBLE PRECISION]
00524         for(jjj=0;jjj<ny2;++jjj)
00525            for(iii=0;iii<nx2;iii++)
00526               Inten[jjj][iii][i] = (double)tempframe[iii+jjj*nx2]; 
00527         drms_free_array(arrin);
00528         }
00529 
00530       free_interpol(&const_param);
00531 
00532       //CREATE OUTPUT ARRAYS
00533       DRMS_Array_t *arrout   = NULL;                                 //array that will contain the phase maps AND Fe I LINE CHARACTERISTICS (LINEWIDTH AND LINEDEPTH)
00534       int axisout[3]         = {4,nx2,ny2};                          //size of the output arrays (nx2=number of columns, ny2=number of rows)
00535       type    = DRMS_TYPE_FLOAT;                                     //type of the output data: FLOAT
00536       arrout  = drms_array_create(type,3,axisout,NULL,&status);      //create the array that will contain the phase maps
00537       if(status != DRMS_SUCCESS || arrout == NULL)
00538       DIE("Error: cannot create array arrout\n");
00539 
00540       DRMS_Array_t *arrout2  = NULL;                                 //array that will contain the reconstructed intensities 
00541       int axisout2[3]        = {nseq,nx2,ny2};                       //size of the output arrays (nx2=number of columns, ny2=number of rows)
00542       type    = DRMS_TYPE_DOUBLE; 
00543       arrout2 = drms_array_create(type,3,axisout2,NULL,&status); 
00544       if(status != DRMS_SUCCESS || arrout2 == NULL)
00545       DIE("Error: cannot create array arrout\n");
00546 
00547       DRMS_Array_t *arrout3  = NULL;                                 //array that will contain the quality flag (did the code converge or not?)
00548       int axisout3[2]        = {nx2,ny2};                            //size of the output arrays (nx2=number of columns, ny2=number of rows)
00549       type    = DRMS_TYPE_SHORT; 
00550       arrout3 = drms_array_create(type,2,axisout3,NULL,&status); 
00551       if(status != DRMS_SUCCESS || arrout3 == NULL)
00552       DIE("Error: cannot create array arrout\n");
00553 
00554       float  *Phig   = arrout->data;                                  //relative phases of the tunable elements
00555       double *IntenR = arrout2->data;                                 //reconstructed intensities
00556       short  *quality= arrout3->data;                                 //quality value
00557       memset(arrout->data,0.0,drms_array_size(arrout));               //initializes the phase values
00558       memset(arrout2->data,0.0,drms_array_size(arrout2));             //initializes the reconstructed intensities values
00559       memset(arrout3->data,0,drms_array_size(arrout3));               //initializes the quality flag values
00560 
00561 
00562 /* MDI Keywords are:
00563 
00564  T_REC      1996.04.29_21:58:16_TAI � 2011.03.25_21:55:44_TAI
00565  ITYPE      "FD", "DARK", ""
00566  T_OBS      1996.04.29_21:58:16_TAI � 2011.03.25_21:55:44_TAI
00567  RECORD     int 0-19
00568  OBSMODE    "STD", "CAL", ""
00569  MDI_CAL1   0, 1, 2, 3, -1
00570  MDI_CAL2   0, 1, 3, -1
00571  MDI_MTM1   range 15 : 114, -1
00572  MDI_MTM2   range 41 : 169, -1
00573  RLC_MTM1   Relative to LC steps or -1
00574  RLC_MTM2   relative to LC steps or -1
00575  TRG_MTM1   e.g. "F1","F2","F3","F4","FC","LC-15",LC",LC+15","FC+1","F1+1","F2+1","F3+1","F4+1"
00576  TRG_MTM2   e.g. "F1","F2","F3","F4","FC","LC-15",LC",LC+15","FC+1","F1+1","F2+1","F3+1","F4+1"
00577 Target and relative correlations:
00578 MTM  Fc F1  F2  F3  F4  LC-15   LC  LC+15
00579 MTM1    -45 -27 -9  9   27  -15 0   15
00580 MTM2    -22 -13 -4  5   14  -15 0   15
00581 NOTE MDI_MTMx = LC(x) + TRG_MTMx
00582  MDI_EXP    1 or 5 or -1
00583  OBSNAME    "cal_detune2b","cal_detune2a","obs_detune2b","obs_detune2a", "obs_detune1b","obs_detune1a"
00584         most are 2a or 2b or missing
00585  MDI_X0
00586  MDI_Y0
00587  CROTA2
00588  CAR_ROT
00589  CRLN_OBS
00590  CRLT_OBS
00591  DSUN_OBS
00592  OBS_VR
00593  OBS_VW
00594  OBS_VN
00595  RSUN_OBS
00596  EARTH_DT
00597  DATAMIN
00598  DATAMAX
00599  DATAMEDN
00600  DATAMEAN
00601  DATARMS
00602  DATASKEW
00603  DATAKURT
00604  DPC        ox404e4400 -> OBSMODE, 0x404d4400 -> CALMODE
00605  T_REC_epoch = 1977.01.01_00:33:12_TAI   ???
00606  T_REC_step = 1
00607 */
00608 
00609 
00610       //READ SOME KEYWORDS OF THE FIRST AND LAST FILTERGRAMS OF THE FULL DETUNE SEQUENCE
00611       Frecnum   = drms_getkey_int(rec2[0],"recnum",&status);          //value of the FSN of the 1st filtergram
00612       if(status != DRMS_SUCCESS)
00613       DIE1("Error: unable to read keyword %s of the first filtergram\n","recnum");
00614 
00615       interntime = drms_getkey_time(rec2[0],"T_OBS",&status);        //value of the T_OBS of the 1st filtergram
00616       if(status != DRMS_SUCCESS)
00617       DIE1("Error: unable to read keyword %s of the first filtergram\n","T_OBS");
00618 
00619       //we take the average value of the radial velocity (in m/s) of observer from 1st and last filtergrams
00620       OBS_VR = drms_getkey_double(rec2[0],VR,&status);            //value of disk center radial velocity for the 1st filtergram
00621       if(status != DRMS_SUCCESS)
00622     {
00623       printf("Error: unable to read the keyword %s of the first filtergram\n",VR);
00624       return(EXIT_FAILURE);
00625     }
00626       OBS_VR += drms_getkey_double(rec2[1],VR,&status);            //value of radial velocity for the last filtergram
00627       if(status != DRMS_SUCCESS)
00628     {
00629       printf("Error: unable to read the keyword %s of the last filtergram\n",VR);
00630       return(EXIT_FAILURE);
00631     }
00632       OBS_VR /= 2.0;
00633 
00634       Lrecnum  = drms_getkey_int(rec2[1],"recnum",&status);          //value of the FSN of the last filtergram
00635       if(status != DRMS_SUCCESS)
00636       DIE1("Error: unable to read keyword %s of the last filtergram\n","recnum");
00637       interntime2= drms_getkey_time(rec2[1],"T_OBS",&status);        //value of the T_OBS of the last filtergram
00638       if(status != DRMS_SUCCESS)
00639       DIE1("Error: unable to read keyword %s of the last filtergram\n","T_OBS");
00640 
00641       printf("recnum value of first filtergram %d\n",Frecnum);          //print value of the FSN of the 1st filtergram
00642       printf("recnum value of last filtergram %d\n",Lrecnum);          //print value of the FSN of the last filtergram
00643 
00644 
00645       //CONVERT THE POSITIONS OF TUNING HCM INTO PHASES (AND CHECK THAT ANGLES ARE 0, 120, OR 240 DEGREES)
00646       for(i=0;i<nseq;++i) 
00647     {
00648       //FOR NB
00649       tuning[i][0]  =  (double)((tuningint[i][0] *   8) % 360);
00650       if(tuning[i][0] != 0.0 && tuning[i][0] != 120.0 && tuning[i][0] != 240.0 && tuning[i][0] != -120.0 && tuning[i][0] != -240.0)
00651         {
00652           printf("Error: detune sequence should have phases of 0, 120, or 240 degrees only: %f\n",tuning[i][0]);
00653           return(EXIT_FAILURE);
00654         }
00655       tuning[i][0]  *=M_PI/180.0;
00656 
00657       //FOR WB
00658       tuning[i][1]  =  (double)((tuningint[i][1] *   8) % 360);
00659       if(tuning[i][1] != 0.0 && tuning[i][1] != 120.0 && tuning[i][1] != 240.0 && tuning[i][1] != -120.0 && tuning[i][1] != -240.0)
00660           {
00661         printf("Error: detune sequence should have phases of 0, 120, or 240 degrees only: %f\n",tuning[i][1]);
00662         return(EXIT_FAILURE);
00663           }
00664       tuning[i][1]  *=M_PI/180.0;
00665 
00666     }
00667 
00668 
00669       //FIND MINIMUM OF RSUN AMONG ALL THE IMAGES OF THE DETUNE SEQUENCE (NB: IN CAL_MODE, RSUN IS NOT THE SOLAR RADIUS)
00670       solarradiusmax = 10000.; //pixels
00671       indexref=0;
00672       for(k=0;k<nseq;k++) if(RSUN[k] < solarradiusmax)
00673     {
00674       solarradiusmax = RSUN[k];
00675       indexref=k;
00676     }
00677       // WARNING PHS - should this be solarradiusmax -= ratio/2 ?? XXXXXXX if problems change back to -= ratio or something...
00678       //solarradiusmax -= ratio/2;  //to avoid problems at the edge, because we average by groups of ratio pixels
00679       solarradiusmax -= ratio;
00680                 // (IF WE DON'T DO THAT, THERE ARE A LOT OF NON-CONVERGENCE AT
00681                                 // THE EDGE OF THE IMAGE DISK, AND THE CODE IS NOTICEABLY SLOWER)
00682       printf("Maximum radius (in pixels) at which phases are computed %f\n",solarradiusmax);
00683 
00684       //MAP OF RADIAL DISTANCES IN PIXELS
00685       printf("Pixel position of the image center used to compute the phase maps: %f %f\n",X0[indexref],Y0[indexref]);
00686       // WARNING PHS - This is where we need to decide if co-registration is required. XXXXXX for obsmode
00687 
00688       for(iii=0;iii<ny*nx;iii++)
00689     {
00690       row    = iii / nx ;
00691       column = iii % nx ;
00692       distance2[iii]=sqrt( ((float)row-Y0[indexref])*((float)row-Y0[indexref])+((float)column-X0[indexref])*((float)column-X0[indexref]) ); 
00693     }
00694 
00695       if (Mode == OBS_MODE)
00696         {
00697         rebinArraySF(distance, distance2, nx2, nx, RSUN[indexref], X0[indexref], Y0[indexref], 512.0);
00698         }
00699       else
00700         {
00701         fresize(&fresizes,distance2,distance,nx,ny,nx,nx2,ny2,nx2,0,0,0.0f);
00702         free_fresize(&fresizes); 
00703         }
00704 
00705 
00706 
00707       /***********************************************************************************************************/
00708       /*BLOCKER FILTER + FRONT WINDOW                                                                            */
00709       /***********************************************************************************************************/
00710  
00711       mdi_filters(nlam, lam, blockerint);
00712 
00713       for(i=0;i<2;++i) FSR[i] = dpi/FSR[i];
00714 
00715 
00716       /***********************************************************************************************************/
00717       /*WAVELENGTH SHIFT DUE TO THE l.o.s. VELOCITY (DEPENDS ON THE LOCATION AT THE SOLAR SURFACE IN OBS_MODE)    */
00718       /* convention: + = away from Sun (redshift)                                                                */
00719       /***********************************************************************************************************/
00720 
00721 
00722       if (Mode == CAL_MODE)
00723         {
00724     for(jjj=0; jjj<ny2; ++jjj)
00725           for(iii=0; iii<nx2; ++iii)
00726         dlam_LOS[jjj*nx2+iii] = dlamdv*OBS_VR;
00727         printf("VELOCITY = %f\n",OBS_VR);
00728         }
00729       else
00730         {
00731         int stat;
00732 // int dims[2] = {nx2,nx2};
00733 // DRMS_Array_t *arr = drms_array_create(DRMS_TYPE_FLOAT,2,dims,NULL,&stat);
00734 // float *arrdata = arr->data;
00735         double DSUN_OBS = drms_getkey_double(rec2[0], "DSUN_OBS", &stat);
00736         double OBS_B0 = drms_getkey_double(rec2[0], "CRLT_OBS", &status); stat += status;
00737         double CROTA2 = drms_getkey_double(rec2[0], "CROTA2", &status); stat += status;
00738         OBS_VW = drms_getkey_double(rec2[0], VW, &status); stat += status;
00739         OBS_VN = (drms_getkey_double(rec2[0], VN, &status)); stat += status;
00740         if (stat || isnan(DSUN_OBS) || isnan(OBS_B0) || isnan(CROTA2) || isnan(OBS_VW) || isnan(OBS_VN))
00741           {
00742       printf("one of CROTA2, CRLT_OBS, DSUN_OBS, OBS_VW, or OBS_VW keywords not available.");
00743       return(EXIT_FAILURE);
00744           }
00745 // printf("OBS_VR=%f, OBS_VW=%f, OBS_VN=%f, CROTA2=%f\n",OBS_VR,OBS_VW,OBS_VN,CROTA2);
00746         double X02 = X0[indexref]/ratio;
00747         double Y02 = Y0[indexref]/ratio;
00748         double rsun2 = RSUN[indexref]/ratio;
00749         double CDELT12 = CDELT1[0]*ratio;
00750         //l.o.s. velocity at pixel location [iii,jjj] 
00751     for(jjj=0; jjj<ny2; ++jjj)
00752           {
00753           for(iii=0; iii<nx2; ++iii)
00754             if (jjj<10 && iii==0)
00755 {
00756           dlam_LOS[jjj*nx2+iii] = dlamdv * OBS_VR; // for bullseye disk averages
00757 // arrdata[jjj*nx2+iii] = 0.0;
00758 }
00759             else
00760           {
00761               double x = iii - X02; // binned pixels
00762               double y = jjj - Y02;
00763           double V_sun = Solar_velocity(x,y,rsun2,OBS_B0,CROTA2,CDELT12,OBS_VR,OBS_VW,OBS_VN);
00764           dlam_LOS[jjj*nx2+iii] = dlamdv * V_sun;
00765 // arrdata[jjj*nx2+iii] = dlamdv * V_sun;
00766           }
00767           }
00768 // DRMS_RecordSet_t *ors = drms_create_records(drms_env, 1, "su_phil.dlam_LOS", DRMS_PERMANENT, &stat);
00769 // DRMS_Record_t *orec = ors->records[0];
00770 // DRMS_Segment_t *oseg = drms_segment_lookupnum(orec,0);
00771 // drms_segment_write(oseg,arr,0);
00772 // drms_close_records(ors,DRMS_INSERT_RECORD);
00773 // return(0);
00774         }
00775 
00776 // SPECIAL BULLS'EYE AVERAGING TEST
00777 // Fill lower left part of rebinned array with averages from 10% to 100% in 10 steps or radius
00778 {
00779 for (jjj=0; jjj<10; jjj++)
00780   dlam_LOS[jjj*nx2+0] = dlamdv*OBS_VR;
00781 
00782 for(i=0;i<nseq;++i) 
00783   {
00784   int ib;
00785   double bullseye[10];
00786   int nbullseye[10];
00787   for (ib=0; ib<10; ib++)
00788     {
00789     bullseye[ib] = 0.0;
00790     nbullseye[ib] = 0;
00791     }
00792   for(jjj=0;jjj<ny2;++jjj)
00793     for(iii=0;iii<nx2;iii++)
00794       {
00795       double val = Inten[jjj][iii][i];
00796       float dist = distance[jjj*nx2 + iii];
00797       if (!isnan(val))
00798         {
00799         for (ib=0; ib<10; ib++)
00800           {
00801           if (dist < round(solarradiusmax*(ib+1.0)/10.0))
00802         {
00803             bullseye[ib] += val;
00804             nbullseye[ib] += 1;
00805             }
00806           }
00807         }
00808       }
00809   for (ib=0; ib<10; ib++)
00810     {
00811     Inten[ib][0][i] = (nbullseye[ib] ? bullseye[ib]/nbullseye[ib] : 0.0);
00812     }
00813   }
00814 }
00815 // END BULLSEYE GENERATION
00816 
00817       
00818       //OPENING OF BINARY FILE CONTAINING CONTRASTS OF TUNABLE ELEMENTS
00819 
00820 XXXXXXX use CON1 and CON2 from mdi_param.c
00821 
00822       
00823       gsl_vector *Residual = NULL;
00824       gsl_matrix *Jac      = NULL;
00825       gsl_matrix *Weights  = NULL;
00826       gsl_matrix *VV       = NULL;
00827       gsl_vector *SS       = NULL;
00828       gsl_vector *work     = NULL;
00829       gsl_vector *tempvec  = NULL;
00830       gsl_vector *tempvec2 = NULL;
00831       
00832       int    location,location1[8];
00833 
00834       /***********************************************************************************************************/
00835       /*LOOP OVER ALL THE PIXELS OF THE CCD                                                                      */
00836       /***********************************************************************************************************/      
00837 
00838       printf("START LOOP\n");
00839       printf("NB: IF THERE ARE MANY NON CONVERGENCE, REMEMBER TO CHECK THE VALUE OF thresh IN HMIparam.h\n");
00840 
00841 #pragma omp parallel default(none) shared(phaseguess,solarradiusmax,factor,Inten,axisout,axisout2,depth0,thresh,width0,nx2,ny2,FSR,CON1,CON2,dpi,lam,tuning,dlam,Phig,distance,blockerint,IntenR,quality,iterstop,shiftw,Icg,fdepthg,fwidthg,a, lam0, dlam_LOS, dlamdv)   private(location,profilef,residual,Residual,i,j,iii,jjj,k,converg,compteur,line,templine,dlinedI0,dlinedw,tempval,profileg,dIntendI0,dIntendw,dIntendIc,dIntendPhi0,dIntendPhi1,dIntendPhi2,tempres,history,err,minimum,Jac,Weights,VV,SS,work,tempvec,tempvec2,jose,phaseE2,phaseE3,phaseE4,phaseE5,contrasts,l,linea,lineb,linec,lined,gaussian,dlinedIc)
00842       {
00843 
00844     Residual = gsl_vector_alloc(nseq);
00845     Jac      = gsl_matrix_alloc(nseq,nparam);
00846     Weights  = gsl_matrix_alloc(nparam,nparam);
00847     VV       = gsl_matrix_alloc(nparam,nparam);
00848     SS       = gsl_vector_alloc(nparam);
00849     work     = gsl_vector_alloc(nparam);
00850     tempvec  = gsl_vector_alloc(nparam);
00851     tempvec2 = gsl_vector_alloc(nparam);
00852     
00853 #pragma omp for
00854     for(jjj=0;jjj<ny2;jjj++) 
00855       {
00856         printf("%d/%d\n",jjj,ny2);
00857         for(iii=0;iii<nx2;iii++) 
00858           {    
00859  
00860         // if(distance[jjj*nx2+iii] <= solarradiusmax) // Changed for BULLSEYE test
00861         if(distance[jjj*nx2+iii] <= solarradiusmax || (jjj<10 && iii==0))
00862           {
00863                   lam0 = dlam_LOS[jjj*nx2+iii];
00864             // for(j=0;j<nseq;j++) printf(" %f ",Inten[jjj][iii][j]);
00865             converg = 0;
00866             compteur= 0;
00867             
00868             Icg[jjj][iii] = thresh;       //estimate of the solar continuum
00869             fdepthg[jjj][iii] = depth0*thresh;//depth of the solar line according to Stenflo & Lindegren (1977)
00870             fwidthg[jjj][iii] = width0;       //width of the solar line according to Stenflo & Lindegren (1977)
00871             
00872             Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]] = (float)((double)phaseguess[0]*M_PI/180./FSR[0]); //guess values for the phases; NB
00873             Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]] = (float)((double)phaseguess[1]*M_PI/180./FSR[1]); // WB
00874             
00875             //NON-TUNABLE HMI TRANSMISSION PROFILE
00876             location    =iii+jjj*nx2;
00877 
00878             //iii = column, jjj = row
00879 
00880 
00881             while(converg == 0)
00882               {
00883             for (k=0;k<nlam;k++)  //BUILD THE SOLAR LINE PROFILE, ASSUMING TWO GAUSSIAN PROFILES
00884               {
00885 
00886                 gaussian[k]= -0.0074*exp(-(lam[k]-lam0+0.200)*(lam[k]-lam0+0.200)/0.13/0.13)-0.021*exp(-(lam[k]-lam0-0.05)*(lam[k]-lam0-0.05)/0.18/0.18); //calibration 13
00887 XXXXXXXXXXX line profile here, shifted for OBS_VR+V
00888 
00889                 dlinedIc[k]= 1.0-gaussian[k];
00890                 l=(lam[k]-lam0)/fwidthg[jjj][iii];
00891                 //if(fabs(lam[k]-lam0) > 1.5) //no need to calculate far from the line center, because it's NaN 
00892                 if(fabs(l) > 26.5) //sinh((x>26.5)^2) = NAN
00893                   {
00894                 line[k]=Icg[jjj][iii]-gaussian[k]*Icg[jjj][iii];
00895                 dlinedw[k] =0.0;
00896                 dlinedI0[k] =0.0;
00897                   }
00898                 else
00899                   {
00900                 line[k]  = Icg[jjj][iii] - fdepthg[jjj][iii] * 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)))-gaussian[k]*Icg[jjj][iii];
00901                 
00902                 //for the derivative relative to fwidthg
00903                 l=(lam[k]-lam0)/(fwidthg[jjj][iii]+0.0001);
00904                 linea[k]= Icg[jjj][iii] - fdepthg[jjj][iii]*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)))-gaussian[k]*Icg[jjj][iii];
00905                 l=(lam[k]-lam0)/(fwidthg[jjj][iii]-0.0001);
00906                 lineb[k]= Icg[jjj][iii]-fdepthg[jjj][iii]*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)))-gaussian[k]*Icg[jjj][iii]; 
00907                 dlinedw[k] = (linea[k]-lineb[k])/(2.0*0.0001);//derivative relative to fwidthg
00908                 
00909                 //for the derivative relative to fdepthg
00910                 l=(lam[k]-lam0)/fwidthg[jjj][iii];
00911                 linec[k]=Icg[jjj][iii]-(fdepthg[jjj][iii]+0.001)*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)))-gaussian[k]*Icg[jjj][iii];
00912                 lined[k]=Icg[jjj][iii]-(fdepthg[jjj][iii]-0.001)*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)))-gaussian[k]*Icg[jjj][iii];
00913                 dlinedI0[k] = (linec[k]-lined[k])/(2.0*0.001);//derivative relative to fdepthg
00914                   }             
00915               }
00916 for (k=0;k<nlam;k++)  //CHECK THE SOLAR LINE PROFILE, etc.
00917 {
00918 int nansfound=0;
00919 if (isnan(blockerint[k])) {nansfound++;printf("nan found in blockerint, k=%d\n",k);}
00920 if (isnan(lam[k])) {nansfound++;printf("nan found in lam, k=%d\n",k);}
00921 if (isnan(linea[k])) {nansfound++;printf("nan found in linea, k=%d\n",k);}
00922 if (isnan(lineb[k])) {nansfound++;printf("nan found in lineb, k=%d\n",k);}
00923 if (isnan(linec[k])) {nansfound++;printf("nan found in linec, k=%d\n",k);}
00924 if (isnan(lined[k])) {nansfound++;printf("nan found in lined, k=%d\n",k);}
00925 if (isnan(profilef[k])) {nansfound++;printf("nan found in profilef, k=%d\n",k);}
00926 if (isnan(gaussian[k])) {nansfound++;printf("nan found in gaussian\n");}
00927 if (nansfound) {printf("jjj=%d, iii=%d, distance=%f\n",jjj,iii,distance[jjj*nx2+iii]);exit(1);}
00928 }
00929 
00930             for(j=0;j<nseq;j++)
00931               {
00932                 residual      = 0.0;
00933                 dIntendI0[j]  = 0.0;
00934                 dIntendw [j]  = 0.0;
00935                 dIntendIc[j]  = 0.0;
00936                 dIntendPhi0[j]= 0.0;
00937                 dIntendPhi1[j]= 0.0;
00938                 
00939                 for (k=0;k<nlam;k++) 
00940                   {
00941                 profileg[k]    = 0.125 * profilef[k] *
00942                                   (1.0+CON1*cos(FSR[0]*(lam[k]+(double)Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][0])) *
00943                   (1.0+CON2*cos(FSR[1]*(lam[k]+(double)Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][1])); 
00944 if (isnan(profileg[k]))
00945 {
00946 printf("nan found in profileg, k=%d, jjj=%d, iii=%d, seq=%d\n",k,jjj,iii,j);
00947 printf("contrasts={%f,%f}, tuning[%d]={%f, %f}\n",contrasts[0],contrasts[1],j,tuning[j][0],tuning[j][1]);
00948 printf("Phig[this jjj,iii]=%f, lam[k]=%f, FSR[0-1]={%f, %f}, phaseguess[0-1]={%f, %f}\n",Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]],lam[k],FSR[0],FSR[1],phaseguess[0],phaseguess[1]);
00949 exit(1);
00950 }
00951 
00952                 residual      += line[k]*profileg[k]*dlam;
00953                 dIntendI0[j]  += (dlinedI0[k]*profileg[k])*dlam;
00954                 dIntendw [j]  += (dlinedw[k] *profileg[k])*dlam;
00955                 dIntendIc[j]  += (dlinedIc[k]*profileg[k])*dlam;
00956                 dIntendPhi0[j]+= (line[k]*(-CON1*FSR[0]*sin(FSR[0]*(lam[k]+(double)Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][0])) * 
00957                     (1.0+CON2*cos(FSR[1]*(lam[k]+(double)Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][1])) / 8.0*profilef[k])*dlam;
00958                 dIntendPhi1[j]+= (line[k]*(-CON2*FSR[1]*sin(FSR[1]*(lam[k]+(double)Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][1])) *
00959                     (1.0+contrasts[0]*cos(FSR[0]*(lam[k]+(double)Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][0])) / 8.0*profilef[k])*dlam;
00960                   }
00961                 jose=Inten[jjj][iii][j]/factor-residual;
00962                 //printf("%d %d %f %f %f %f %f\n",iii,jjj,dIntendI0[j],dIntendIc[j],dIntendw[j],dIntendPhi0[j],dIntendPhi1[j]);
00963                 gsl_vector_set(Residual,j,jose);                               //residual is computed     
00964                 gsl_matrix_set(Jac,(size_t)j,0,dIntendI0[j]);                  //Jacobian matrix
00965                 gsl_matrix_set(Jac,(size_t)j,1,dIntendIc[j]);
00966                 gsl_matrix_set(Jac,(size_t)j,2,dIntendw[j]);
00967                 gsl_matrix_set(Jac,(size_t)j,3,dIntendPhi0[j]);
00968                 gsl_matrix_set(Jac,(size_t)j,4,dIntendPhi1[j]);
00969                 gsl_matrix_set(Jac,(size_t)j,5,dIntendPhi2[j]);
00970                 
00971               }//for(j=0;j<nseq;++j)
00972 
00973 // printf("D\n");
00974             //SVD algorithm from the GNU Scientific Library
00975             //A is MxN. A = U S V^T for M >= N. On output the matrix A is replaced by U.
00976             gsl_linalg_SV_decomp(Jac,VV,SS,work);
00977             gsl_matrix_set_zero(Weights);
00978             for(i=0;i<nparam;i++) gsl_matrix_set(Weights,(size_t)i,(size_t)i,1.0/gsl_vector_get(SS,(size_t)i)); //creates the diagonal matrix
00979             gsl_blas_dgemv(CblasTrans  ,1.0,Jac,Residual   ,0.0,tempvec );
00980             gsl_blas_dgemv(CblasNoTrans,1.0,Weights,tempvec,0.0,tempvec2);
00981             gsl_blas_dgemv(CblasNoTrans,1.0,VV,tempvec2    ,0.0,tempvec );
00982             
00983             //compute relative changes in the parameters
00984             tempres[0]   = fabs(gsl_vector_get(tempvec,0)/fdepthg[jjj][iii]);
00985             tempres[1]   = fabs(gsl_vector_get(tempvec,1)/Icg[jjj][iii]);
00986             tempres[2]   = fabs(gsl_vector_get(tempvec,2)/fwidthg[jjj][iii]);
00987             tempres[3]   = fabs(gsl_vector_get(tempvec,3)/(double)Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]);
00988             tempres[4]   = fabs(gsl_vector_get(tempvec,4)/(double)Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]);
00989             
00990 
00991 // printf("D0: tempres[x], x0: %g, X1: %g, X2: %g, X3: %g, X4: %g\n",tempres[0],tempres[1],tempres[2],tempres[3],tempres[4]);
00992             err[compteur]= fmax(tempres[0],fmax(tempres[1],fmax(tempres[2],fmax(tempres[3],fmax(tempres[4])))));
00993                           
00994 // printf("D1: compteur=%d, err=%g, iterstop=%g\n",compteur,err[compteur],iterstop);
00995             if( (compteur == maxsteps-1) && (err[compteur] > iterstop) ) converg = 2; //no convergence 
00996             if(err[compteur] <= iterstop)                                converg = 1; //convergence (we stop the iteration when the maximum relative change is less than iterstop)
00997 // printf("E, converg=%d\n",converg);
00998             
00999             if(converg == 2)
01000               {
01001                 gsl_sort_smallest_index(&minimum,1,err,1,maxsteps);                //finds the minimum of err
01002                 fdepthg[jjj][iii] = history[minimum][0];
01003                 Icg[jjj][iii]     = history[minimum][1];
01004                 fwidthg[jjj][iii] = history[minimum][2];
01005                 Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]] = (float)history[minimum][3];
01006                 Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]] = (float)history[minimum][4];
01007                 quality[location]=1;                                               //the quality flag is raised to mention that the code did not converge
01008               }
01009             
01010             if(converg != 2)
01011               {
01012                 fdepthg[jjj][iii] += gsl_vector_get(tempvec,0);
01013                 if( (fdepthg[jjj][iii] < (0.2*thresh*depth0)) || (fdepthg[jjj][iii] > (3.0*depth0*thresh)) || isnan(fdepthg[jjj][iii])) fdepthg[jjj][iii] = depth0*thresh;
01014                 Icg[jjj][iii]     += gsl_vector_get(tempvec,1);
01015                 if( (Icg[jjj][iii] < (0.2*thresh)) || (Icg[jjj][iii] > (3.0*thresh)) || isnan(Icg[jjj][iii])) Icg[jjj][iii] = thresh;
01016                 fwidthg[jjj][iii] += gsl_vector_get(tempvec,2);
01017                 if( (fwidthg[jjj][iii] > (5.0*width0)) || (fwidthg[jjj][iii] < (0.3*width0)) || isnan(fwidthg[jjj][iii])) fwidthg[jjj][iii] = width0;
01018                 Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]] += (float)gsl_vector_get(tempvec,3);
01019                 if( fabsf(Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]) > (1.2*M_PI/FSR[0]) ) Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]] = 0.;
01020                 Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]] += (float)gsl_vector_get(tempvec,4);
01021                 if( fabsf(Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]) > (1.2*M_PI/FSR[1]) ) Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]] = 0.;
01022               }
01023 // printf("G, converg=%d\n",converg);
01024             
01025             history[compteur][0]=fdepthg[jjj][iii];
01026             history[compteur][1]=Icg[jjj][iii];
01027             history[compteur][2]=fwidthg[jjj][iii];
01028             history[compteur][3]=(double)Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]];
01029             history[compteur][4]=(double)Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]];
01030                        
01031             compteur += 1;
01032             
01033             if(converg == 2)
01034                           {
01035               // for (j=maxsteps-4; j<maxsteps;j++)
01036               j = maxsteps-1;
01037                           printf("no convergence %d,%d: %d: fd=%f, Ic=%f, fw=%f, NB=%f, WB=%f, E1=%f\n",jjj,iii,j,history[j][0],history[j][1],history[j][2],history[j][3]*FSR[0]*Deg2Rad,history[j][4]*FSR[1]*Deg2Rad,history[j][5]*FSR[2]*Deg2Rad);
01038                           }
01039               }//while(converg==0)
01040         
01041 
01042 
01043             /***********************************************************************************************************/
01044             /*RECONSTRUCTING THE DETUNE SEQUENCE                                                                       */
01045             /***********************************************************************************************************/       
01046   
01047             for(k=0;k<nlam;++k)
01048               {
01049             l=(lam[k]-lam0)/fwidthg[jjj][iii];
01050             //if(fabs(lam[k]-lam0) <= 1.5) //no need to calculate far from the line center
01051             if(fabs(l) <= 26.5)
01052               {
01053                 line[k]  = Icg[jjj][iii]-fdepthg[jjj][iii]*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)))-gaussian[k]*Icg[jjj][iii];                
01054               }
01055             else line[k]=Icg[jjj][iii]-gaussian[k]*Icg[jjj][iii];
01056               }
01057             for(j=0;j<nseq;j++)
01058               {
01059             residual   = 0.0;
01060             for(k=0;k<nlam;++k)
01061               {
01062                 profileg[k]=profilef[k]*0.125 *
01063                               (1.0+contrasts[0]*cos(FSR[0]*(lam[k]+(double)Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][0])) *
01064                               (1.0+contrasts[1]*cos(FSR[1]*(lam[k]+(double)Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]])+tuning[j][1]));
01065                 residual      += line[k]*profileg[k]*dlam;
01066               }
01067             IntenR[j+iii*axisout2[0]+jjj*axisout2[0]*axisout2[1]] = residual;
01068               }
01069             
01070             Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=(float)((double)Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]*(FSR[0]*180.0/M_PI));
01071             Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=(float)((double)Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]*(FSR[1]*180.0/M_PI));
01072             Phig[2+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=fwidthg[jjj][iii];//we save sigma in Angstroms
01073             Phig[3+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=fdepthg[jjj][iii]/Icg[jjj][iii];//we save the linedepth, for a continuum of 1
01074           }//if(distance <= solarradiusmax)
01075         else
01076           {
01077             Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0; //NB Michelson
01078             Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0; //WB Michelson
01079             Icg[jjj][iii]=0.0;
01080             fwidthg[jjj][iii]=0.0;
01081             fdepthg[jjj][iii]=0.0;
01082             Phig[2+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0;
01083             Phig[3+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0;
01084           }
01085           }//for jjj
01086       }//for iii
01087     
01088     
01089     gsl_vector_free(Residual);
01090     gsl_matrix_free(Jac);
01091     gsl_matrix_free(VV);
01092     gsl_vector_free(SS);
01093     gsl_vector_free(work);
01094     gsl_vector_free(tempvec);
01095     gsl_vector_free(tempvec2);
01096     gsl_matrix_free(Weights);     
01097       }//pragma omp parallel    
01098       
01099 
01100       /***********************************************************************************************************/
01101       /*COMPUTE THE SPATIAL AVERAGES OF THE DETUNE SEQUENCES (ORIGINAL AND RECONSTRUCTED)                        */
01102       /***********************************************************************************************************/       
01103       
01104       
01105       double meanInten[nseq], meanIntenR[nseq];
01106       double NBphase=0.0, WBphase=0.0, E1phase=0.0;
01107       double meanIcg=0.0,meanwidthg=0.0,meandepthg=0.0;
01108       compteur = 0;
01109 
01110       for(j=0;j<nseq;++j)
01111     {
01112       meanInten[j] =0.0;
01113       meanIntenR[j]=0.0;
01114     }
01115 
01116       for(jjj=0;jjj<ny2;jjj++) 
01117     {
01118       for(iii=0;iii<nx2;iii++) 
01119         {       
01120           //if(distance[jjj*nx2+iii] <= (solarradiustable/CDELT1[indexref]) && quality[iii+jjj*nx2] == 0)
01121           if(distance[jjj*nx2+iii] <= (900.0/CDELT1[indexref]) && quality[iii+jjj*nx2] == 0) //900 arcsec instead of the full disk 976"
01122         {
01123           NBphase += Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]];
01124           WBphase += Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]];
01125           meanIcg += Icg[jjj][iii];
01126           meanwidthg += fwidthg[jjj][iii]; 
01127           meandepthg += fdepthg[jjj][iii]; 
01128           compteur+= 1;
01129           for(j=0;j<nseq;++j)
01130             {
01131               // meanInten[j] +=Inten[jjj][iii][j]/factor;                             //original intensities
01132               meanInten[j] +=Inten[jjj][iii][j]/nseq;                                  //original intensities
01133               meanIntenR[j]+=factor/nseq*IntenR[j+iii*axisout2[0]+jjj*axisout2[0]*axisout2[1]];    //reconstructed intensities
01134             }
01135         }
01136         }
01137     }
01138       
01139       NBphase = NBphase/(double)compteur;
01140       WBphase = WBphase/(double)compteur;
01141       meanIcg = meanIcg/(double)compteur;
01142       meanwidthg = meanwidthg/(double)compteur;
01143       meandepthg = meandepthg/(double)compteur;
01144 
01145 
01146       drms_free_array(arrout2);
01147 
01148 
01149 
01150       /***********************************************************************************************************/
01151       /*ATTEMPT TO CORRECT THE NON-CONVERGENCE                                                                   */
01152       /***********************************************************************************************************/      
01153 
01154       printf("CORRECTION OF PIXELS WHERE THERE WAS NO CONVERGENCE, USING NEIGHBOR AVERAGING\n");
01155       for(jjj=0;jjj<ny2;jjj++) 
01156     {
01157       for(iii=0;iii<nx2;iii++) 
01158         {
01159           location=iii+jjj*nx2;
01160           if(quality[location] ==1) //there was no convergence
01161         {
01162           //location of 8 neighboring pixels
01163           location1[0]=(jjj+1)*nx2+iii;     //(jjj+1,iii) (row,column)
01164           location1[1]=(jjj+1)*nx2+(iii+1); //(jjj+1,iii+1)
01165           location1[2]=jjj*nx2+(iii+1);     //(jjj,iii+1)
01166           location1[3]=(jjj-1)*nx2+(iii+1); //(jjj-1,iii+1)
01167           location1[4]=(jjj-1)*nx2+iii;     //(jjj-1,iii)
01168           location1[5]=(jjj-1)*nx2+(iii-1); //(jjj-1,iii-1)
01169           location1[6]=jjj*nx2+(iii-1);     //(jjj,iii-1)
01170           location1[7]=(jjj+1)*nx2+(iii-1); //(jjj+1,iii-1)
01171           compteur=0;
01172           Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0;
01173           Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0;
01174           Phig[2+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0;
01175           Phig[3+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=0.0;
01176           for(j=0;j<8;++j)
01177             {
01178               if(quality[location1[j]] == 0) //there was convergence
01179             {
01180               Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]+=Phig[0+axisout[0]*(location1[j])];
01181               Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]+=Phig[1+axisout[0]*(location1[j])];
01182               Phig[2+iii*axisout[0]+jjj*axisout[0]*axisout[1]]+=Phig[2+axisout[0]*(location1[j])];
01183               Phig[3+iii*axisout[0]+jjj*axisout[0]*axisout[1]]+=Phig[3+axisout[0]*(location1[j])];
01184               compteur+=1;
01185             }         
01186             }
01187           if(compteur != 0)
01188             {
01189               Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=Phig[0+iii*axisout[0]+jjj*axisout[0]*axisout[1]]/(float)compteur;
01190               Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=Phig[1+iii*axisout[0]+jjj*axisout[0]*axisout[1]]/(float)compteur;
01191               Phig[2+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=Phig[2+iii*axisout[0]+jjj*axisout[0]*axisout[1]]/(float)compteur;
01192               Phig[3+iii*axisout[0]+jjj*axisout[0]*axisout[1]]=Phig[3+iii*axisout[0]+jjj*axisout[0]*axisout[1]]/(float)compteur;
01193             }
01194         }
01195         }
01196     }
01197 
01198       /***********************************************************************************************************/
01199       /*DISPLAY RESULTS                                                                                          */
01200       /***********************************************************************************************************/      
01201  
01202       printf("RESULTS\n");
01203       printf("-----------------------------------------------------------------------------\n");
01204       printf("\n SPATIALLY AVERAGED INTENSITIES \n");
01205       for(j=0;j<nseq;++j) printf("%f %f\n",meanInten[j],meanIntenR[j]);
01206       printf("\n SPATIALLY AVERAGED PHASES (IN DEGREES)\n");
01207       printf("NB Michelson: %f\n",NBphase);
01208       printf("WB Michelson: %f\n",WBphase);
01209       printf("\n SPATIALLY AVERAGED WIDTH, DEPTH, AND CONTINUUM\n");
01210       printf("WIDTH: %f \n",meanwidthg);
01211       printf("DEPTH: %f \n",meandepthg);
01212       printf("CONTINUUM: %f \n",meanIcg);
01213       printf("\n COTUNE TABLE \n");
01214       printf("WB    NB\n");
01215       int    table[2][20];
01216       
01217       cotunetable(NBphase,WBphase,table);
01218       
01219       for(j = 0; j < 20; ++j)
01220       {
01221         printf("%d   %d\n",table[0][j],table[1][j]);      
01222       } 
01223 
01224       /***********************************************************************************************************/
01225       /*CREATE A RECORD FOR THE PHASE MAPS                                                                       */
01226       /***********************************************************************************************************/      
01227       
01228       DRMS_RecordSet_t *dataout = NULL;
01229       DRMS_Record_t  *recout    = NULL;
01230       DRMS_Segment_t *segout    = NULL;
01231 
01232       dataout = drms_create_records(drms_env,1,dsout,DRMS_PERMANENT,&status);
01233       
01234       if (status != DRMS_SUCCESS)
01235     {
01236       printf("Could not create a record for the phase maps\n");
01237       return(EXIT_FAILURE);
01238     }
01239       if (status == DRMS_SUCCESS)
01240     {
01241     printf("Writing a record on the DRMS for the phase maps\n");
01242     recout = dataout->records[0];
01243 
01244       //WRITE KEYWORDS
01245       
01246 XXXXXXXXXXXXX fix keywords
01247     status += drms_setkey_int(recout,"START_rec",Frecnum);              //FSN_START
01248     status += drms_setkey_time(recout,"T_START",interntime);          //T_START
01249     status += drms_setkey_int(recout,"STOP_rec",Lrecnum); 
01250     status += drms_setkey_time(recout,"T_STOP",interntime2);          //T_STOP
01251     Mrecnum= (Frecnum+Lrecnum)/2;                           //FSN_REC (PRIME KEY)
01252     interntime3 = (interntime+interntime2)/2.0;                  //T_REC   (PRIME KEY)
01253     status += drms_setkey_int(recout,"RECNUM",keyvalue3);
01254     status += drms_setkey_time(recout,"T_REC",interntime3);
01255     status += drms_setkey_string(recout,COMMENTS,COMMENT);        //COMMENT ON THE DATA 
01256     status += drms_setkey_float(recout,CRPIX1,X0[indexref]+1.0);  //CRPIX1
01257     status += drms_setkey_float(recout,CRPIX2,Y0[indexref]+1.0);  //CRPIX2
01258     status += drms_setkey_float(recout,SCALE,CDELT1[indexref]);   //CDELT1
01259     status += drms_setkey_int(recout,FOCUS,HCFTID[indexref]);     //HCFTID
01260     status += drms_setkey_float(recout,RADIUS,solarradiusmax);    //R_SUN (NOT NECESSARILY SOLAR RADIUS, BUT SIZE OF THE IMAGE. solarradiusmax differs from RSUN[indexref] by -ratio pixels)
01261     status += drms_setkey_int(recout,HCMNB0,(int)round(-NBphase/6.0)+60); // NB motor position for average phase
01262     status += drms_setkey_int(recout,HCMWB0,(int)round(-WBphase/6.0)+60); // WB motor position for average phase
01263     status += drms_setkey_int(recout,NXs,nx2);
01264     char  DATEOBS[256];
01265     sprint_time(DATEOBS,CURRENT_SYSTEM_TIME,"UTC",1);
01266     status += drms_setkey_string(recout,DATEs,DATEOBS);           //DATE AT WHICH THE FILE WAS CREATED
01267     status += drms_setkey_float(recout,FSRNBs,2.0*M_PI/FSR[0]);
01268     status += drms_setkey_float(recout,FSRWBs,2.0*M_PI/FSR[1]);
01269     status += drms_setkey_float(recout,CBLOCKERs,centerblocker2);
01270         status += drms_setkey_float(recout, "PHASENBM", NBphase);
01271         status += drms_setkey_float(recout, "PHASEWBM", WBphase);
01272         status += drms_setkey_float(recout, "LW_MEAN", meanwidthg);
01273         status += drms_setkey_float(recout, "LD_MEAN", meandepthg);
01274         status += drms_setkey_float(recout, "CON_MEAN", meanIcg);
01275         status += drms_setkey_string(recout,"MODE",(Mode==OBS_MODE ? "OBS_MODE" : "CAL_MODE"));
01276     if (status)
01277           {
01278           printf("Failed to write one or more keywords.");
01279           return(EXIT_FAILURE);
01280           }
01281 
01282       //WRITE DATA SEGMENTS
01283       segout = drms_segment_lookupnum(recout, 0);
01284       drms_segment_write(segout, arrout, 0);                       //write the file containing the phase maps
01285       segout = drms_segment_lookupnum(recout, 1);
01286       drms_segment_write(segout, arrout3,0);                       //write the file containing the quality flag
01287           printf("Phases and quality segments written\n");
01288           // Write phasemaps as separate files
01289           char *segnames[] = {"NB", "WB", "LW", "LD"};
01290           for (k=0; k<4; k++)
01291             {
01292 // printf("Segment %s to be written.\n",segnames[k]);
01293             int status;
01294             int axis[2]         = {ny2,nx2};           //size of the output arrays (nx2=number of columns, ny2=number of rows)
01295             DRMS_Array_t *arr = drms_array_create(DRMS_TYPE_FLOAT,2,axis,NULL,&status);  //create the array that will contain the phase maps
01296             float *data = arr->data;
01297             for (jjj=0; jjj<ny2; jjj++)
01298               for (iii=0; iii<nx2; iii++)
01299         data[jjj*nx2 + iii] = Phig[k+iii*axisout[0]+jjj*axisout[0]*axisout[1]];
01300         segout = drms_segment_lookup(recout, segnames[k]);
01301 if (!segout){printf("lookup for %s segment failed\n", segnames[k]);continue;}
01302         status = drms_segment_write(segout, arr,0);                       //write the file containing a phasemap or value
01303         printf("Segment for %s %s written.\n",segnames[k],(status? "FAILED TO BE" : ""));
01304             drms_free_array(arr);
01305             }
01306 
01307       //CLOSE THE RECORDS
01308       drms_close_records(dataout, DRMS_INSERT_RECORD);             //insert the record in DRMS
01309     }      
01310       drms_free_array(arrout);   
01311       drms_free_array(arrout3); 
01312     }
01313     
01314   return error;
01315   
01316 }
01317 
01318 
01319 // Functions needed for OBS_MODE LOS velocity
01320 
01321 double Solar_velocity(double x,double y,double rsun2,double OBS_B0,double CROTA2, double CDELT, double obs_vr, double obs_vw, double obs_vn)
01322   {
01323   double Msun = 200.0;   // Mean velocity
01324   double A = 1995.0;
01325   double B = -315.0;
01326   double LS0 = -455.0;
01327   double LS1 = 140.0;
01328   double LS2 = 720.0;
01329   OBS_B0 *= Deg2Rad;
01330   CROTA2 *= Deg2Rad;
01331   double arcsec2Rad = Deg2Rad/3600;
01332   double P = -CROTA2;
01333   double cosa = cos(CROTA2);
01334   double sina = sin(CROTA2);
01335   double X = x*cosa - y*sina; // pixels
01336   double Y = y*cosa + x*sina;
01337   double r = sqrt(x*x + y*y)/rsun2;
01338   double lon, sinlat, coslat, sig;
01339   double sinB2, rfac, r2, z, val;
01340   if (r < 1.0)
01341       img2sphere(X/rsun2, Y/rsun2, rsun2*CDELT*arcsec2Rad, OBS_B0, 0.0, P, NULL, NULL, &lon, &sinlat, &coslat, &sig, NULL, NULL);
01342   else
01343     return(0.0);
01344   // away from Sun is positive, away from East is positive, away from South is positive, red shift is positive.
01345   double obs_v = obs_vr*cos(sig) - obs_vw*sin(X*CDELT*arcsec2Rad) - obs_vn*sin(Y*CDELT*arcsec2Rad);
01346   sinB2 = sinlat * sinlat;
01347   rfac = coslat * cos(OBS_B0) * sin(lon);
01348   r2 = r * r;
01349   z = 1 - sqrt(1.0 - r2); // 1 - costheta
01350 //printf("x=%f, y=%f, lon=%f, lat=%f\n",X/rsun2,Y/rsun2,lon/Deg2Rad,asin(sinlat)/Deg2Rad);
01351   // val = Msun + rfac * (A + B* sinB2 * (1.0 + sinB2)) + z * (LS0 + z * (LS1 + z * LS2)) ;
01352   // val = obs_vr;
01353   // val = obs_v;
01354   // val = obs_vr + Msun + rfac * A;
01355   val = obs_v + Msun + rfac * (A + B* sinB2 * (1.0 + sinB2)) + z * (LS0 + z * (LS1 + z * LS2)) ;
01356   return(val);
01357   }
01358 
01359 /*
01360  *                     FUNCTIONS
01361  *
01362  *          img2sphere()
01363  *
01364  *  Written: Rick Bogart
01365  *
01366  */
01367 
01368 int img2sphere (double x, double y, double ang_r, double latc, double lonc,
01369     double pa, double *p_rho, double *p_lat, double *p_lon, double *p_sinlat,
01370     double *p_coslat, double *p_sig, double *p_mu, double *p_chi)
01371   {
01372   /*
01373    *  Taken from cartography.c, written by Rick Bogart
01374    *  Map projected coordinates (x, y) to (lon, lat) and (rho | sig, chi)
01375    *
01376    *  Arguments:
01377    *    x            Plate locations, in units of the image radius, relative
01378    *    y                to the image center
01379    *                      These x,y are 0:1 not arcsecs and not rotated by crota2.
01380    *    ang_r         Apparent semi-diameter of sun (angular radius of sun at
01381    *                      the observer's tangent line)
01382    *    latc          Latitude of disc center, uncorrected for light travel time
01383    *    lonc          Longitude of disc center
01384    *    pa            Position angle of solar north on image, measured eastward
01385    *                      from north (sky coordinates)
01386    *  Return values:
01387    *    rho           Angle point:sun_center:observer
01388    *    lon           Heliographic longitude
01389    *    lat           Heliographic latitude
01390    *    sinlat        sine of heliographic latitude
01391    *    coslat        cosine of heliographic latitude
01392    *    sig           Angle point:observer:sun_center
01393    *    mu            cosine of angle between the point:observer line and the
01394    *                      local normal
01395    *    chi           Position angle on image measured westward from solar
01396    *                      north
01397    *
01398    *  All angles are in radians.
01399    *  Return value is 1 if point is outside solar radius (in which case the
01400    *    heliographic coordinates and mu are meaningless), 0 otherwise.
01401    *  It is assumed that the image is direct; the x or y coordinates require a
01402    *    sign change if the image is inverted.
01403    *
01404    * modified to ignore setting return values for NULL pointers.
01405    *
01406    */
01407 
01408   double w_rho, w_lat, w_lon, w_sinlat, w_coslat, w_sig, w_mu, w_chi;
01409   double *rho= &w_rho, *lat= &w_lat, *lon= &w_lon, *sinlat= &w_sinlat, *coslat= &w_coslat, *sig= &w_sig, *mu= &w_mu, *chi= &w_chi;
01410   static double ang_r0 = 0.0, sinang_r = 0.0, tanang_r = 0.0;
01411   static double latc0 = 0.0, coslatc = 1.0, sinlatc = 0.0;
01412   double cosr, sinr, sinlon, sinsig;
01413 
01414   if (p_rho) rho=p_rho;
01415   if (p_lat) lat=p_lat;
01416   if (p_lon) lon=p_lon;
01417   if (p_sinlat) sinlat=p_sinlat;
01418   if (p_coslat) coslat=p_coslat;
01419   if (p_sig) sig=p_sig;
01420   if (p_mu) mu=p_mu;
01421   if (p_chi) chi=p_chi;
01422 
01423   if (ang_r != ang_r0) {
01424       sinang_r = sin(ang_r);
01425       tanang_r = tan(ang_r);
01426       ang_r0 = ang_r;
01427     }
01428 
01429   if (latc != latc0) {
01430       sinlatc = sin(latc);
01431       coslatc = cos(latc);
01432       latc0 = latc;
01433     }
01434 
01435   *chi = atan2 (x, y) + pa;
01436   while (*chi > 2 * M_PI) *chi -= 2 * M_PI;
01437   while (*chi < 0.0) *chi += 2 * M_PI;
01438   /*  Camera curvature correction, no small angle approximations  */
01439   *sig = atan (hypot (x, y) * tanang_r);
01440   sinsig = sin (*sig);
01441   *rho = asin (sinsig / sinang_r) - *sig;
01442   if (*sig > ang_r) return (-1);
01443   *mu = cos (*rho + *sig);
01444   sinr = sin (*rho);
01445   cosr = cos (*rho);
01446   *sinlat = sinlatc * cosr + coslatc * sinr * cos (*chi);
01447   *coslat = sqrt (1.0 - *sinlat * *sinlat);
01448   *lat = asin (*sinlat);
01449   sinlon = sinr * sin (*chi) / *coslat;
01450   *lon = asin (sinlon);
01451   if (cosr < (*sinlat * sinlatc)) *lon = M_PI - *lon;
01452   *lon += lonc;
01453   while (*lon < 0.0) *lon += 2 * M_PI;
01454   while (*lon >= 2 * M_PI) *lon -= 2 * M_PI;
01455   return (0);
01456   }
01457 
01458 int rebinArraySF(float *outData, float *inData, int outNx, int inNx, float rsun, float x0, float y0, float fillvalue)
01459   { //  Rebin only using data inside the Sun, else use fillvalue
01460   float rsun2 = 0.999*rsun*rsun; // require whole pixels for 4K images
01461   int nmiss = 0;
01462   float *inp, *outp;
01463   int inNy = inNx;
01464   int outNy = outNx;
01465   int inI, inJ, outI, outJ, i;
01466   int binsize = inNx/outNx;
01467   double inRow[inNx], outRow[outNx];
01468   int inRowN[inNx], outRowN[outNx];
01469 
01470   for (outJ=0; outJ < outNy; outJ++)
01471     {
01472     int inRow0 = outJ * binsize;
01473     for (outI=0; outI<outNx; outI++)
01474       {
01475       outRowN[outI] = 0;
01476       outRow[outI] = 0.0;
01477       }
01478     for (inJ = inRow0; inJ < inRow0+binsize; inJ++)
01479       {
01480       inp = inData + inJ*inNx;
01481       for (outI=0; outI<outNx; outI++)
01482         for (i=0; i<binsize; i++, inp++)
01483           {
01484           inI = outI*binsize + i;
01485           if (*inp != DRMS_MISSING_FLOAT && (((inI-x0)*(inI-x0) + (inJ-y0)*(inJ-y0)) < rsun2))
01486             {
01487             outRow[outI] += *inp;
01488             outRowN[outI]++;
01489             }
01490           }
01491        }
01492     for (outI=0; outI < outNx; outI++)
01493       if (outRowN[outI])
01494         outData[outJ*outNx + outI] = outRow[outI]/outRowN[outI];
01495       else
01496         {
01497         outData[outJ*outNx + outI] = fillvalue;
01498         nmiss++;
01499         }
01500     }  // end outJ
01501   }
proj/lev1.5_hmi/apps/MDI_phasemaps.c
