JSOC_Documentation: proj/lev1.5_hmi/libs/lev15/Dopplergram2.c Source File

00001 /*-----------------------------------------------------------------------------------------*/
00002 /*                                                                                         */
00003 /* Program to compute Dopplergrams using a MDI-like algorithm                              */
00004 /*(compute the phases of the 1st and 2nd Fourier coefficients)                             */
00005 /* Author: S. Couvidat (based on a code by J. Schou)                                       */
00006 /* Version 1.8 April 20, 2009                                                              */
00007 /*                                                                                         */
00008 /* USES BICUBIC INTERPOLATION UNLIKE Dopplergram.c                                         */
00009 /* uses a MDI-like algorithm with 5 or 6 tuning positions                                  */
00010 /* averages the velocities returned by 1st and 2nd Fourier                                 */
00011 /* coefficient, and by LCP and RCP                                                         */
00012 /* the code also estimates the Fe I linewidth, linedepth, and                              */
00013 /* the continuum intensity                                                                 */
00014 /*                                                                                         */
00015 /* FOR THE MEANING OF ROWS AND COLUMNS WE FOLLOW THE CONVENTION:                           */
00016 /* A 2D ARRAY A IN C IS DEFINED AS A[row][column]                                          */
00017 /*                                                                                         */
00018 /* CONVENTION FOR THE VELOCITIES:                                                          */
00019 /* POSITIVE VELOCITIES CORRESPOND TO REDSHIFT (MOVEMENTS AWAY FROM THE OBSERVER)           */
00020 /*-----------------------------------------------------------------------------------------*/
00021 
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <math.h>
00025 #include <complex.h>
00026 #include <jsoc_main.h>
00027 #include <omp.h>                      //OpenMP header
00028 
00029 #undef I                              //I is the complex number (0,1). We un-define it to avoid confusion
00030 
00031 struct parameterDoppler {             //structure to provide some parameters defined in HMIparam.h to Dopplergram()
00032   double FSRNB;
00033   double FSRWB;
00034   double FSRE1;
00035   double FSRE2;
00036   double FSRE3;
00037   double FSRE4;
00038   double FSRE5;
00039   double dlamdv;
00040   int maxVtest;
00041   int maxNx;
00042   int ntest;
00043   double dvtest;
00044   float MISSINGDATA;
00045   float MISSINGRESULT;
00046 };
00047 
00048 
00049 /*-------------------------------------------------------------------*/
00050 /* bicubic interpolation routine derived from Numerical Recipes in C */
00051 /*-------------------------------------------------------------------*/
00052 
00053 void bcucof(float y[4],float y1[4],float y2[4],float y12[4],float d1,float d2,float c[4][4])
00054  {
00055    static int wt[16][16]= {
00056      1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
00057      0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
00058      -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
00059      2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
00060      0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
00061      0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
00062      0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
00063      0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
00064      -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
00065      0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
00066      9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
00067      -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
00068      2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
00069      0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
00070      -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
00071      4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1};
00072    
00073    int l,k,j,i;
00074    float xx,d1d2,cl[16],x[16];
00075 
00076    d1d2=d1*d2;
00077    for(i=0;i<=3;i++){
00078      x[i]=y[i];
00079      x[i+4]=y1[i]*d1;
00080      x[i+8]=y2[i]*d2;
00081      x[i+12]=y12[i]*d1d2;
00082    }
00083    for(i=0;i<=15;i++){
00084      xx=0.0;
00085      for(k=0;k<=15;k++) xx += wt[i][k]*x[k];
00086      cl[i]=xx;
00087    }
00088    l=0;
00089    for(i=0;i<=3;i++)
00090      for(j=0;j<=3;j++) c[i][j]=cl[l++];
00091 
00092  }
00093 
00094 
00095 /*-----------------------------------------------------------------------------------------*/
00096 /* MDI-like algorithm                                                                      */
00097 /*-----------------------------------------------------------------------------------------*/
00098 
00099 
00100 int Dopplergram2(DRMS_Array_t **arrLev1p,DRMS_Array_t **arrLev15,int framelistSize,DRMS_Array_t *Lookuptable,float Rsun,float X0,float Y0,struct parameterDoppler DopplerParameters,int *MISSVALS,int *SATVALS,float cdelt1)
00101 {
00102 
00103   double FSR[7];
00104   double dlamdv;
00105   int maxVtest;
00106   int maxNx;
00107   int ntest;
00108   double dvtest;
00109   float MISSINGDATA;
00110   float MISSINGRESULT;
00111   int  i,j,iii,k;                                                      //loop variables
00112   float ExtraCrop = 50.;//we crop the data at Rsun+ExtraCrop in pixels
00113 
00114   double magnetic = 1.0/(2.0*4.67E-5*0.000061733433*2.5*299792458.0); //Lande factor=2.5 for Fe I line
00115 
00116   //parameters of the Fe I line profiles observed by Roger Ulrich (OBTAINED BY CODE centertolimb.pro; MAKE SURE IT'S THE SAME AS IN lookup.c)
00117   double cost[3]={1.0,0.70710678,0.5};                    //cos(theta) where theta is the angular distance from disk center
00118   double minimumCoeffs[2]={0.41922611,0.24190794};        //minimum intensity, assuming the continuum is at 1 (result from a Gaussian fit in the range [-0.055,0.055] Angstroms)
00119   double FWHMCoeffs[2]={151.34559,-58.521771};            //FWHM of the solar line in Angstrom (result from a Gaussian fit in the range [-0.055,0.055] Angstroms)
00120   double FWHM,minimum,angulardistance,minimumR,minimumL;
00121 
00122   FSR[0]=DopplerParameters.FSRNB;
00123   FSR[1]=DopplerParameters.FSRWB;
00124   FSR[2]=DopplerParameters.FSRE1;
00125   FSR[3]=DopplerParameters.FSRE2;
00126   FSR[4]=DopplerParameters.FSRE3;
00127   FSR[5]=DopplerParameters.FSRE4;
00128   FSR[6]=DopplerParameters.FSRE5;
00129   dlamdv=DopplerParameters.dlamdv;
00130   maxVtest=DopplerParameters.maxVtest;
00131   maxNx=DopplerParameters.maxNx;
00132   ntest=DopplerParameters.ntest;
00133   dvtest=DopplerParameters.dvtest;
00134   MISSINGDATA=DopplerParameters.MISSINGDATA;
00135   MISSINGRESULT=DopplerParameters.MISSINGRESULT;
00136   double vtest[ntest];                                                //internally calculations are done in double precision
00137   float poly[ntest],poly2[ntest];                                     //for the spatial interpolation of the look-up tables (THE TABLES ARE ASSUMED TO BE OF TYPE FLOAT)
00138 
00139   int  status      = 0; 
00140   int  error       = 0;                                               //error code returned by the routine: 0=no error, !=0 means error
00141   int  nRows,nColumns;
00142   DRMS_Type_t type;
00143   float distance;
00144 
00145   //check whether or not framelistSize is an even number
00146   if(framelistSize != 10 && framelistSize != 12 && framelistSize != 16 && framelistSize != 20)
00147     {
00148       printf("Error in subroutine computing the Dopplergrams: the framelist size is not 10, 12, 16, or 20\n");
00149       error=3;
00150       return error;
00151     }
00152 
00153   int N=framelistSize/2;                                             //number of different wavelengths
00154   double *cosi;
00155   double *sini;
00156   double *cos2i;
00157   double *sin2i;
00158   cosi = (double *)malloc(N*sizeof(double));  
00159   sini = (double *)malloc(N*sizeof(double));  
00160   cos2i= (double *)malloc(N*sizeof(double));  
00161   sin2i= (double *)malloc(N*sizeof(double));  
00162   if(cosi == NULL || sini == NULL || cos2i == NULL || sin2i == NULL)
00163     {
00164       printf("Error: memory could not be allocated in the Dopplergram() function\n");
00165       exit(EXIT_FAILURE);
00166     }
00167   float L[N],R[N]; //NO MALLOCS BECAUSE THEY ARE WRITTEN INSIDE THE OMP LOOP, THE LEVEL 1P DATA ARE ASSUMED TO BE OF TYPE FLOAT
00168  
00169 
00170   //check that all level 1p data have the same data type (IS THAT NECESSARY?)
00171   for(i=0;i<framelistSize;++i)
00172     {
00173       type    = arrLev1p[i]->type;                                  //float if produced by Jesper's routine
00174       if (type != DRMS_TYPE_FLOAT)                                  //if the type is not FLOAT
00175     {
00176       printf("Error in subroutine computing the Dopplergrams: data type of level 1p data is not FLOAT\n");
00177       error = 1;
00178       free(cosi);
00179       free(sini);
00180       free(cos2i);
00181       free(sin2i);
00182       return error;
00183     }
00184     }
00185   
00186   //with C convention Array[row][column], axis[0] is the number of columns 
00187   nRows    = arrLev1p[0]->axis[1];
00188   nColumns = arrLev1p[0]->axis[0];
00189 
00190   //check that all filtergrams have the same dimensions (IS THAT NECESSARY?)
00191   for(i=1;i<framelistSize;++i)
00192     {
00193       if (arrLev1p[i]->axis[1] != nRows || arrLev1p[i]->axis[0] != nColumns )//if filtergrams are not nRows*nColumns
00194     {
00195       printf("Error in subroutine computing the Dopplergrams: dimensions of level 1p data are not %d x %d \n",nRows,nColumns);
00196       error = 2;
00197       free(cosi);
00198       free(sini);
00199       free(cos2i);
00200       free(sin2i);
00201       return error;
00202     }
00203     }
00204 
00205 
00206   //THE LEVEL 1.5 DATA ARE ASSUMED TO BE OF TYPE FLOAT
00207   float *lam0g  = arrLev15[0]->data ;                                 //Dopplergram
00208   float *B0g    = arrLev15[1]->data;                                  //magnetogram
00209   float *Idg    = arrLev15[2]->data;                                  //linedepth
00210   float *widthg = arrLev15[3]->data;                                  //linewidth
00211   float *I0g    = arrLev15[4]->data;                                  //continuum
00212 
00213   memset(lam0g,  0.0, drms_array_size(arrLev15[0]));                  //fill the observable arrays with 0
00214   memset(B0g  ,  0.0, drms_array_size(arrLev15[1]));
00215   memset(Idg  ,  0.0, drms_array_size(arrLev15[2]));
00216   memset(widthg, 0.0, drms_array_size(arrLev15[3]));
00217   memset(I0g   , 0.0, drms_array_size(arrLev15[4]));
00218   
00219   
00220   //variables for the MDI-like algorithm
00221   double FSRNB   = FSR[0];                       //FSR Narrow-Band Michelson, in Angstroms (WILL CHANGE ONCE THE VALUE IS ACCURATELY MEASURED)
00222   double dtune   = FSRNB/2.5;                    //wavelength separation between each tuning position, nominally 68.8 mA  (SHOULD BE THE SAME FOR 5 OR 6 WAVELENGTHS)
00223   double dv      = 1.0/dlamdv;                   //conversion factor from wavelength to velocity
00224   double dvtune  = dtune*dv;
00225   double *tune   = NULL, *angle=NULL;
00226   tune=(double *)malloc(N*sizeof(double));
00227   if(tune == NULL)
00228     {
00229       printf("Error: unable to allocate memory to tune\n");
00230       exit(EXIT_FAILURE);
00231     }
00232   angle=(double *)malloc(N*sizeof(double));
00233   if(angle == NULL)
00234     {
00235       printf("Error: unable to allocate memory to angle\n");
00236       exit(EXIT_FAILURE);
00237     }
00238 
00239 
00240   if(N == 6)
00241     {
00242       printf("USING 6 WAVELENGTHS\n");
00243       tune[0]=+2.5;//I0 (location of the filter)
00244       tune[1]=+1.5;//I1
00245       tune[2]=+0.5;//I2
00246       tune[3]=-0.5;//I3
00247       tune[4]=-1.5;//I4
00248       tune[5]=-2.5;//I5
00249       for(i=0;i<N;++i) angle[i]=tune[i];
00250     }
00251   if(N == 5)
00252     {
00253       printf("USING 5 WAVELENGTHS\n");
00254       tune[0]=+2.0;
00255       tune[1]=+1.0;
00256       tune[2]= 0.0;
00257       tune[3]=-1.0;
00258       tune[4]=-2.0;
00259       for(i=0;i<N;++i) angle[i]=tune[i];
00260     }
00261   if(N == 8)
00262     {
00263       printf("USING 8 WAVELENGTHS\n");
00264       tune[7]=+3.5;//I7
00265       tune[0]=+2.5;//I0 (location of the filter)
00266       tune[1]=+1.5;//I1
00267       tune[2]=+0.5;//I2
00268       tune[3]=-0.5;//I3
00269       tune[4]=-1.5;//I4
00270       tune[5]=-2.5;//I5
00271       tune[6]=-3.5;//I6
00272       for(i=0;i<N;++i) angle[i]=tune[i];
00273     }
00274   if(N == 10)
00275     {
00276       printf("USING 10 WAVELENGTHS\n");
00277       tune[9]=+4.5;//I9
00278       tune[7]=+3.5;//I7
00279       tune[0]=+2.5;//I0 (location of the filter)
00280       tune[1]=+1.5;//I1
00281       tune[2]=+0.5;//I2
00282       tune[3]=-0.5;//I3
00283       tune[4]=-1.5;//I4
00284       tune[5]=-2.5;//I5
00285       tune[6]=-3.5;//I6
00286       tune[8]=-4.5;//I7
00287       for(i=0;i<N;++i) angle[i]=tune[i];
00288     }
00289 
00290   double period = (double)(N-1)*dtune;
00291   double pv1,pv2;
00292   double f1LCPc,f1RCPc,f1LCPs,f1RCPs,f2LCPc,f2RCPc,f2LCPs,f2RCPs;
00293   double vLCP,vRCP,v2LCP,v2RCP;
00294   double temp,temp2,temp3,tempbis,temp2bis,temp3bis;
00295   double meanL=0.0,meanR=0.0;
00296 
00297   pv1 = dvtune*(double)(N-1);
00298   pv2 = pv1/2.;
00299   
00300   for(i=0;i<N;++i) 
00301     {
00302       tune[i] = tune[i]*dtune;
00303       angle[i]= angle[i]*2.0*M_PI/(double)N;
00304       cosi[i] = cos(angle[i]);
00305       sini[i] = sin(angle[i]);
00306       cos2i[i]= cos(2.0*angle[i]);
00307       sin2i[i]= sin(2.0*angle[i]);
00308     }
00309   
00310   
00311   //array containing the look-up table (type FLOAT)
00312   float *lookupt = Lookuptable->data ;
00313   int axist[3];                  //dimensions of the look-up tables
00314   axist[0]=Lookuptable->axis[0]; //number of input velocities (ntest*2)
00315   axist[1]=Lookuptable->axis[1]; //number of columns (256)    
00316   axist[2]=Lookuptable->axis[2]; //number of rows (256)
00317   printf("Dimensions of the look-up tables: %d %d %d\n",axist[0],axist[1],axist[2]);
00318 
00319   if(axist[0] > maxVtest || axist[1] > maxNx || axist[2] > maxNx)
00320     {
00321       printf("Error in subroutine computing the Dopplergrams: dimensions of the look-up tables exceed what is allowed: %d %d %d\n",maxVtest,maxNx,maxNx); //if there is a problem
00322       error = 4;
00323       free(cosi);
00324       free(sini);
00325       free(cos2i);
00326       free(sin2i);
00327       free(tune);
00328       free(angle);
00329       return error;        
00330     }
00331 
00332   //to accommodate the bilinear interpolation code we format the 1D lookup tables into 2d arrays
00333   float lookupt2D[Lookuptable->axis[1]][Lookuptable->axis[2]];
00334 
00335 
00336   //we reconstruct the test (input) velocities used to produce the look-up tables
00337   //WARNING: MUST BE THE SAME AS IN lookup.c
00338   for(i=0;i<ntest;++i) vtest[i] = dvtest*((double)i-((double)ntest-1.0)/2.0); 
00339 
00340   int   index_lo=0;
00341   int   index_hi=ntest-1;
00342   float RR1,RR2;     //variables for the bilinear interpolation of the look-up tables
00343   int   x0,y0,x1,y1;
00344   float xa,xb,ya,yb;
00345   int   ratio,indexL,indexR,indexL2,indexR2,row,column,step=10;
00346   long  loc0,loc1,loc2,loc3,loc4,loc5,loc6,loc7,loc8,loc9,loc10,loc11,loc12,loc13,loc14,loc15; //coded 8 bytes
00347   double Kfourier;
00348   float *tempvec=NULL;
00349   Kfourier = dtune/period*2.0;
00350   ratio = nRows/axist[2];
00351   float minlookupt1=0.0,maxlookupt1=0.0,minlookupt2=0.0,maxlookupt2=0.0;
00352   
00353   int MISSVALS2=0; //the number of missing values (NaN)
00354   int SATVALS2 =0; //number of staurated values
00355   float coeffs[4][4];
00356   float u,t,yy[4],yy1[4],yy2[4],yy12[4];
00357   float offset=((float)ratio-1.0)/2.0; //because the phase maps were rebinned using rebin() which means that pixel 0 is actually (ratio-1)/2 on the initial grid
00358 
00359   /***********************************************************************************************************/
00360   /*LOOP OVER ALL THE PIXELS OF THE FILTERGRAMS                                                              */
00361   /***********************************************************************************************************/
00362 
00363 #pragma omp parallel default(none) reduction(+:MISSVALS2,SATVALS2) shared(step,arrLev1p,cosi,sini,cos2i,sin2i,pv1,pv2,index_lo,index_hi,vtest,period,dtune,dv,I0g,B0g,Idg,lam0g,widthg,magnetic,axist,ratio,lookupt,nRows,nColumns,MISSINGDATA,MISSINGRESULT,Kfourier,Rsun,X0,Y0,ntest,tune,N,cost,minimumCoeffs,FWHMCoeffs,offset,ExtraCrop,cdelt1) private(tempvec,iii,L,R,f1LCPc,f1RCPc,f1LCPs,f1RCPs,vLCP,vRCP,f2LCPc,f2RCPc,f2LCPs,f2RCPs,temp,tempbis,temp2,temp2bis,temp3,temp3bis,meanL,meanR,v2LCP,v2RCP,x0,y0,x1,y1,RR1,RR2,i,loc0,loc1,loc2,loc3,loc4,xa,xb,ya,yb,indexL,indexR,indexL2,indexR2,poly,poly2,row,column,distance,j,minlookupt1,maxlookupt1,minlookupt2,maxlookupt2,FWHM,minimum,angulardistance,minimumR,minimumL,coeffs,u,t,k,yy,yy1,yy2,yy12,loc5,loc6,loc7,loc8,loc9,loc10,loc11,loc12,loc13,loc14,loc15)
00364  {
00365 
00366 #pragma omp for
00367    for(iii=0;iii<nRows*nColumns;++iii)
00368     {
00369       //with the convention adopted for a 2D array, the index iii is defined as iii=column+row*nColumns
00370       row   =iii / nColumns;
00371       column=iii % nColumns;
00372 
00373       distance = sqrt(((float)row-Y0)*((float)row-Y0)+((float)column-X0)*((float)column-X0)); //distance in pixels
00374       if(distance <= (Rsun+ExtraCrop))
00375         {
00376 
00377           //CALCULATE THE SOLAR LINE PARAMETER AT THE DISTANCE FROM DISK CENTER (USING A LAW DERIVED FROM THE 3 Fe I PROFILES PROVIDED BY ROGER ULRICH)
00378           //if(distance <= Rsun) angulardistance=cos(asin(distance/Rsun)); //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)
00379           //else angulardistance=0.0;
00380           //FWHM=FWHMCoeffs[0]+FWHMCoeffs[1]*angulardistance;
00381           //FWHM=FWHM/2.0/sqrt(log(2.0))/1000.; //we convert from FWHM in mA to sigma in A
00382           //minimum=minimumCoeffs[0]+minimumCoeffs[1]*angulardistance; //minimum intensity Id
00383 
00384           //CALCULATE THE SOLAR LINE PARAMETER AT THE DISTANCE FROM DISK CENTER (USING A LAW DERIVED FROM A LINEWIDTH MAP OBTAINED WITH THE MDI-LIKE ALGORITHM)
00385           distance=distance*cdelt1; //convert from pixels to arcsecs
00386           FWHM=100.67102+0.015037016*distance-0.00010128197*distance*distance+3.1548385E-7*distance*distance*distance-3.7298102E-10*distance*distance*distance*distance+1.7275788E-13*distance*distance*distance*distance*distance;
00387           FWHM=FWHM/2.0/sqrt(log(2.0))/1000.; //we convert from FWHM in mA to sigma in A
00388 
00389           /*-------------------------------------------------------------*/
00390           /* MDI-like algorithm                                          */
00391           /*-------------------------------------------------------------*/
00392 
00393           for(i=0;i<N;++i) 
00394             {
00395               tempvec=(float *)arrLev1p[i*2]->data;   //LCP (ASSUMES THE LEVEL 1p DATA ARE STORED LCP/RCP BACK-TO-BACK) IN THE ORDER I0, I1, I2, I3, I4, and I5
00396               L[i]= tempvec[iii] ;
00397               tempvec=(float *)arrLev1p[i*2+1]->data; //RCP
00398               R[i]= tempvec[iii] ;            
00399             }
00400 
00401           //First and Second Fourier coefficients
00402           f1LCPc=0.0;
00403           f1RCPc=0.0;
00404           f1LCPs=0.0;
00405           f1RCPs=0.0;
00406           f2LCPc=0.0;
00407           f2RCPc=0.0;
00408           f2LCPs=0.0;
00409           f2RCPs=0.0;
00410           for(i=0;i<N;++i) 
00411             {
00412               f1LCPc += cosi[i] *(double)L[i];
00413               f1RCPc += cosi[i] *(double)R[i];
00414               f1LCPs += sini[i] *(double)L[i];
00415               f1RCPs += sini[i] *(double)R[i];
00416               f2LCPc += cos2i[i]*(double)L[i];
00417               f2RCPc += cos2i[i]*(double)R[i];
00418               f2LCPs += sin2i[i]*(double)L[i];
00419               f2RCPs += sin2i[i]*(double)R[i];
00420             }
00421           
00422           vLCP    = atan2(-f1LCPs,-f1LCPc)*pv1/2.0/M_PI; //-f1LCPs and -f1LCPc so that the jump is at 180 degrees and not 0 degrees
00423           vRCP    = atan2(-f1RCPs,-f1RCPc)*pv1/2.0/M_PI;
00424               v2LCP   = atan2(-f2LCPs,-f2LCPc)*pv2/2.0/M_PI;
00425           v2RCP   = atan2(-f2RCPs,-f2RCPc)*pv2/2.0/M_PI;
00426           v2LCP   = fmod((v2LCP-vLCP+10.5*pv2),pv2)-pv2/2.0+vLCP; //we use the uncorrected velocity, i.e. phase, of the 1st Fourier coefficient to correct for the estimate of v2LCP and v2RCP, because the range of velocities obtained with the second Fourier coefficient is half the range of the first Fourier coefficient
00427           v2RCP   = fmod((v2RCP-vRCP+10.5*pv2),pv2)-pv2/2.0+vRCP; 
00428 
00429           if(isnan(f1LCPc) || isnan(f1RCPc) || isnan(f1LCPs) || isnan(f1RCPs))
00430             {
00431               MISSVALS2+=1;
00432               lam0g[iii]  = MISSINGRESULT;
00433               B0g[iii]    = MISSINGRESULT;
00434               widthg[iii] = MISSINGRESULT;
00435               Idg[iii]    = MISSINGRESULT;
00436               I0g[iii]    = MISSINGRESULT;
00437               continue;
00438             }
00439 
00440           /*-------------------------------------------------------------*/
00441           /* bicubic interpolation of the look-up tables at pixel (x,y)  */
00442           /*-------------------------------------------------------------*/
00443                   
00444           //NB: it depends on how the filtergrams rebinning from 4096*4096 to axist[1]*axist[2] was done in phasemaps.c
00445           //find the 4 neighbors (x0,y0), (x0,y1), (x1,y0), and (x1,y1) of (column,row) on the grid of the look-up tables, and deal with boundary problems
00446 
00447           x0 = floor( ((float)column-offset)/(float)ratio); //a column number
00448           y0 = floor( ((float)row-offset)/(float)ratio);    //a row number
00449           x1 = x0+1;
00450           y1 = y0+1;
00451 
00452           //boundary issues
00453           if(x1 >= axist[1]-1)
00454             {
00455               x0 = x0-2;
00456               x1 = x1-2;
00457             }
00458           if(y1 >= axist[2]-1)
00459             {
00460               y0 = y0-2;
00461               y1 = y1-2;
00462             }
00463           if(x1 < 2)
00464             {
00465               x0 = 1;
00466               x1 = 2;
00467             }
00468           if(y1 < 2)
00469             {
00470               y0 = 1;
00471               y1 = 2;
00472             }
00473 
00474           xb = (((float)column-offset)-(float)x0*(float)ratio)/(float)ratio;
00475           yb = (((float)row-offset)-(float)y0*(float)ratio)/(float)ratio;
00476 
00477           //for bicubic interpolation
00478           loc0= x0   *axist[0]+ y0   *axist[0]*axist[1]; //x0
00479           loc1= x1   *axist[0]+ y0   *axist[0]*axist[1]; //x1
00480           loc2= x1   *axist[0]+ y1   *axist[0]*axist[1]; //x2
00481           loc3= x0   *axist[0]+ y1   *axist[0]*axist[1]; //x3
00482           loc4=(x0-1)*axist[0]+ y0   *axist[0]*axist[1]; //x4
00483           loc5= x0   *axist[0]+(y0-1)*axist[0]*axist[1]; //x5
00484           loc6= x1   *axist[0]+(y0-1)*axist[0]*axist[1]; //x6
00485           loc7=(x1+1)*axist[0]+ y0   *axist[0]*axist[1]; //x7
00486           loc8=(x1+1)*axist[0]+ y1   *axist[0]*axist[1]; //x8
00487           loc9= x1   *axist[0]+(y1+1)*axist[0]*axist[1]; //x9
00488           loc10= x0   *axist[0]+(y1+1)*axist[0]*axist[1]; //x10
00489           loc11=(x0-1)*axist[0]+ y1   *axist[0]*axist[1]; //x11
00490           loc12=(x0-1)*axist[0]+(y0-1)*axist[0]*axist[1]; //x12
00491           loc13=(x1+1)*axist[0]+(y0-1)*axist[0]*axist[1]; //x13
00492           loc14=(x1+1)*axist[0]+(y1+1)*axist[0]*axist[1]; //x14
00493           loc15=(x0-1)*axist[0]+(y1+1)*axist[0]*axist[1]; //x15
00494 
00495           indexL =ntest+1;
00496           indexR =ntest+1;
00497           indexL2=ntest+1;
00498           indexR2=ntest+1;
00499 
00500           //initializes arrays (NOT NEEDED ACTUALLY)
00501           memset(poly,0.0,sizeof(poly));
00502           memset(poly2,0.0,sizeof(poly2));
00503 
00504 
00505           //minimum velocity
00506           i=0;
00507           yy[0]=lookupt[loc0+i];
00508           yy[1]=lookupt[loc1+i];
00509           yy[2]=lookupt[loc2+i];
00510           yy[3]=lookupt[loc3+i];
00511           yy1[0]=(lookupt[loc1+i]-lookupt[loc4+i])/2.0; //derivative in x direction
00512           yy1[1]=(lookupt[loc7+i]-lookupt[loc0+i])/2.0;
00513           yy1[2]=(lookupt[loc8+i]-lookupt[loc3+i])/2.0;
00514           yy1[3]=(lookupt[loc2+i]-lookupt[loc11+i])/2.0;
00515           yy2[0]=(lookupt[loc3+i] -lookupt[loc5+i])/2.0; //derivative in y direction
00516           yy2[1]=(lookupt[loc2+i] -lookupt[loc6+i])/2.0;
00517           yy2[2]=(lookupt[loc9+i] -lookupt[loc1+i])/2.0;
00518           yy2[3]=(lookupt[loc10+i]-lookupt[loc0+i])/2.0;
00519           yy12[0]=(lookupt[loc2+i] -lookupt[loc6+i] -lookupt[loc11+i]+lookupt[loc12+i])/4.0; //cross-derivatives
00520           yy12[1]=(lookupt[loc8+i] -lookupt[loc13+i]-lookupt[loc3+i] +lookupt[loc5+i] )/4.0;
00521           yy12[2]=(lookupt[loc14+i]-lookupt[loc7+i] -lookupt[loc10+i]+lookupt[loc0+i] )/4.0;
00522           yy12[3]=(lookupt[loc9+i] -lookupt[loc1+i] -lookupt[loc15+i]+lookupt[loc4+i] )/4.0;
00523           bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00524           u=yb;
00525           t=xb;
00526           poly[0]=0.0;
00527           for (i = 3; i >= 0; i--) poly[0] = t*poly[0] + ((coeffs[i][3]*u+coeffs[i][2])*u + coeffs[i][1])*u + coeffs[i][0];
00528 
00529           i=ntest;
00530           yy[0]=lookupt[loc0+i];
00531           yy[1]=lookupt[loc1+i];
00532           yy[2]=lookupt[loc2+i];
00533           yy[3]=lookupt[loc3+i];
00534           yy1[0]=(lookupt[loc1+i]-lookupt[loc4+i])/2.0; //derivative in x direction
00535           yy1[1]=(lookupt[loc7+i]-lookupt[loc0+i])/2.0;
00536           yy1[2]=(lookupt[loc8+i]-lookupt[loc3+i])/2.0;
00537           yy1[3]=(lookupt[loc2+i]-lookupt[loc11+i])/2.0;
00538           yy2[0]=(lookupt[loc3+i] -lookupt[loc5+i])/2.0; //derivative in y direction
00539           yy2[1]=(lookupt[loc2+i] -lookupt[loc6+i])/2.0;
00540           yy2[2]=(lookupt[loc9+i] -lookupt[loc1+i])/2.0;
00541           yy2[3]=(lookupt[loc10+i]-lookupt[loc0+i])/2.0;
00542           yy12[0]=(lookupt[loc2+i] -lookupt[loc6+i] -lookupt[loc11+i]+lookupt[loc12+i])/4.0; //cross-derivatives
00543           yy12[1]=(lookupt[loc8+i] -lookupt[loc13+i]-lookupt[loc3+i] +lookupt[loc5+i] )/4.0;
00544           yy12[2]=(lookupt[loc14+i]-lookupt[loc7+i] -lookupt[loc10+i]+lookupt[loc0+i] )/4.0;
00545           yy12[3]=(lookupt[loc9+i] -lookupt[loc1+i] -lookupt[loc15+i]+lookupt[loc4+i] )/4.0;
00546           bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00547           u=yb;
00548           t=xb;
00549           poly2[0]=0.0;
00550           for (i = 3; i >= 0; i--) poly2[0] = t*poly2[0] + ((coeffs[i][3]*u+coeffs[i][2])*u + coeffs[i][1])*u + coeffs[i][0];
00551 
00552           minlookupt1 = poly[0];
00553           minlookupt2 = poly2[0];
00554 
00555 
00556           //maximum velocity
00557           i=ntest-1;
00558           yy[0]=lookupt[loc0+i];
00559           yy[1]=lookupt[loc1+i];
00560           yy[2]=lookupt[loc2+i];
00561           yy[3]=lookupt[loc3+i];
00562           yy1[0]=(lookupt[loc1+i]-lookupt[loc4+i])/2.0; //derivative in x direction
00563           yy1[1]=(lookupt[loc7+i]-lookupt[loc0+i])/2.0;
00564           yy1[2]=(lookupt[loc8+i]-lookupt[loc3+i])/2.0;
00565           yy1[3]=(lookupt[loc2+i]-lookupt[loc11+i])/2.0;
00566           yy2[0]=(lookupt[loc3+i] -lookupt[loc5+i])/2.0; //derivative in y direction
00567           yy2[1]=(lookupt[loc2+i] -lookupt[loc6+i])/2.0;
00568           yy2[2]=(lookupt[loc9+i] -lookupt[loc1+i])/2.0;
00569           yy2[3]=(lookupt[loc10+i]-lookupt[loc0+i])/2.0;
00570           yy12[0]=(lookupt[loc2+i] -lookupt[loc6+i] -lookupt[loc11+i]+lookupt[loc12+i])/4.0; //cross-derivatives
00571           yy12[1]=(lookupt[loc8+i] -lookupt[loc13+i]-lookupt[loc3+i] +lookupt[loc5+i] )/4.0;
00572           yy12[2]=(lookupt[loc14+i]-lookupt[loc7+i] -lookupt[loc10+i]+lookupt[loc0+i] )/4.0;
00573           yy12[3]=(lookupt[loc9+i] -lookupt[loc1+i] -lookupt[loc15+i]+lookupt[loc4+i] )/4.0;
00574           bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00575           u=yb;
00576           t=xb;
00577           poly[ntest-1]=0.0;
00578           for (i = 3; i >= 0; i--) poly[ntest-1] = t*poly[ntest-1] + ((coeffs[i][3]*u+coeffs[i][2])*u + coeffs[i][1])*u + coeffs[i][0];
00579 
00580           i=2*ntest-1;
00581           yy[0]=lookupt[loc0+i];
00582           yy[1]=lookupt[loc1+i];
00583           yy[2]=lookupt[loc2+i];
00584           yy[3]=lookupt[loc3+i];
00585           yy1[0]=(lookupt[loc1+i]-lookupt[loc4+i])/2.0; //derivative in x direction
00586           yy1[1]=(lookupt[loc7+i]-lookupt[loc0+i])/2.0;
00587           yy1[2]=(lookupt[loc8+i]-lookupt[loc3+i])/2.0;
00588           yy1[3]=(lookupt[loc2+i]-lookupt[loc11+i])/2.0;
00589           yy2[0]=(lookupt[loc3+i] -lookupt[loc5+i])/2.0; //derivative in y direction
00590           yy2[1]=(lookupt[loc2+i] -lookupt[loc6+i])/2.0;
00591           yy2[2]=(lookupt[loc9+i] -lookupt[loc1+i])/2.0;
00592           yy2[3]=(lookupt[loc10+i]-lookupt[loc0+i])/2.0;
00593           yy12[0]=(lookupt[loc2+i] -lookupt[loc6+i] -lookupt[loc11+i]+lookupt[loc12+i])/4.0; //cross-derivatives
00594           yy12[1]=(lookupt[loc8+i] -lookupt[loc13+i]-lookupt[loc3+i] +lookupt[loc5+i] )/4.0;
00595           yy12[2]=(lookupt[loc14+i]-lookupt[loc7+i] -lookupt[loc10+i]+lookupt[loc0+i] )/4.0;
00596           yy12[3]=(lookupt[loc9+i] -lookupt[loc1+i] -lookupt[loc15+i]+lookupt[loc4+i] )/4.0;
00597           bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00598           u=yb;
00599           t=xb;
00600           poly2[ntest-1]=0.0;
00601           for (i = 3; i >= 0; i--) poly2[ntest-1] = t*poly2[ntest-1] + ((coeffs[i][3]*u+coeffs[i][2])*u + coeffs[i][1])*u + coeffs[i][0];
00602 
00603           maxlookupt1 = poly[ntest-1];
00604           maxlookupt2 = poly2[ntest-1];
00605 
00606           for(i=step;i<ntest;i=i+step) //make sure (ntest-1)/step is an integer
00607             {
00608               //poly[i]  =ya*(lookupt[loc1+i]*xa+lookupt[loc2+i]*xb)+yb*(lookupt[loc3+i]*xa+lookupt[loc4+i]*xb);                          //for 1st Fourier coefficient
00609               //poly2[i] =ya*(lookupt[loc1+i+ntest]*xa+lookupt[loc2+i+ntest]*xb)+yb*(lookupt[loc3+i+ntest]*xa+lookupt[loc4+i+ntest]*xb);  //for 2nd Fourier coefficient
00610 
00611               yy[0]=lookupt[loc0+i];
00612               yy[1]=lookupt[loc1+i];
00613               yy[2]=lookupt[loc2+i];
00614               yy[3]=lookupt[loc3+i];
00615               yy1[0]=(lookupt[loc1+i]-lookupt[loc4+i])/2.0; //derivative in x direction
00616               yy1[1]=(lookupt[loc7+i]-lookupt[loc0+i])/2.0;
00617               yy1[2]=(lookupt[loc8+i]-lookupt[loc3+i])/2.0;
00618               yy1[3]=(lookupt[loc2+i]-lookupt[loc11+i])/2.0;
00619               yy2[0]=(lookupt[loc3+i] -lookupt[loc5+i])/2.0; //derivative in y direction
00620               yy2[1]=(lookupt[loc2+i] -lookupt[loc6+i])/2.0;
00621               yy2[2]=(lookupt[loc9+i] -lookupt[loc1+i])/2.0;
00622               yy2[3]=(lookupt[loc10+i]-lookupt[loc0+i])/2.0;
00623               yy12[0]=(lookupt[loc2+i] -lookupt[loc6+i] -lookupt[loc11+i]+lookupt[loc12+i])/4.0; //cross-derivatives
00624               yy12[1]=(lookupt[loc8+i] -lookupt[loc13+i]-lookupt[loc3+i] +lookupt[loc5+i] )/4.0;
00625               yy12[2]=(lookupt[loc14+i]-lookupt[loc7+i] -lookupt[loc10+i]+lookupt[loc0+i] )/4.0;
00626               yy12[3]=(lookupt[loc9+i] -lookupt[loc1+i] -lookupt[loc15+i]+lookupt[loc4+i] )/4.0;
00627               bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00628               u=yb;
00629               t=xb;
00630               poly[i]=0.0;
00631               for (j = 3; j >= 0; j--) poly[i] = t*poly[i] + ((coeffs[j][3]*u+coeffs[j][2])*u + coeffs[j][1])*u + coeffs[j][0];
00632 
00633               yy[0]=lookupt[loc0+i+ntest];
00634               yy[1]=lookupt[loc1+i+ntest];
00635               yy[2]=lookupt[loc2+i+ntest];
00636               yy[3]=lookupt[loc3+i+ntest];
00637               yy1[0]=(lookupt[loc1+i+ntest]-lookupt[loc4+i+ntest])/2.0; //derivative in x direction
00638               yy1[1]=(lookupt[loc7+i+ntest]-lookupt[loc0+i+ntest])/2.0;
00639               yy1[2]=(lookupt[loc8+i+ntest]-lookupt[loc3+i+ntest])/2.0;
00640               yy1[3]=(lookupt[loc2+i+ntest]-lookupt[loc11+i+ntest])/2.0;
00641               yy2[0]=(lookupt[loc3+i+ntest] -lookupt[loc5+i+ntest])/2.0; //derivative in y direction
00642               yy2[1]=(lookupt[loc2+i+ntest] -lookupt[loc6+i+ntest])/2.0;
00643               yy2[2]=(lookupt[loc9+i+ntest] -lookupt[loc1+i+ntest])/2.0;
00644               yy2[3]=(lookupt[loc10+i+ntest]-lookupt[loc0+i+ntest])/2.0;
00645               yy12[0]=(lookupt[loc2+i+ntest] -lookupt[loc6+i+ntest] -lookupt[loc11+i+ntest]+lookupt[loc12+i+ntest])/4.0; //cross-derivatives
00646               yy12[1]=(lookupt[loc8+i+ntest] -lookupt[loc13+i+ntest]-lookupt[loc3+i+ntest] +lookupt[loc5+i+ntest] )/4.0;
00647               yy12[2]=(lookupt[loc14+i+ntest]-lookupt[loc7+i+ntest] -lookupt[loc10+i+ntest]+lookupt[loc0+i+ntest] )/4.0;
00648               yy12[3]=(lookupt[loc9+i+ntest] -lookupt[loc1+i+ntest] -lookupt[loc15+i+ntest]+lookupt[loc4+i+ntest] )/4.0;
00649               bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00650               u=yb;
00651               t=xb;
00652               poly2[i]=0.0;
00653               for (j = 3; j >= 0; j--) poly2[i] = t*poly2[i] + ((coeffs[j][3]*u+coeffs[j][2])*u + coeffs[j][1])*u + coeffs[j][0];
00654 
00655               if(poly[i] > vLCP && poly[i-step]   <= vLCP)
00656             {
00657               for(j=i-step+1;j<=i;j++)
00658                 {
00659 
00660                   yy[0]=lookupt[loc0+j];
00661                   yy[1]=lookupt[loc1+j];
00662                   yy[2]=lookupt[loc2+j];
00663                   yy[3]=lookupt[loc3+j];
00664                   yy1[0]=(lookupt[loc1+j]-lookupt[loc4+j])/2.0; //derivative in x direction
00665                   yy1[1]=(lookupt[loc7+j]-lookupt[loc0+j])/2.0;
00666                   yy1[2]=(lookupt[loc8+j]-lookupt[loc3+j])/2.0;
00667                   yy1[3]=(lookupt[loc2+j]-lookupt[loc11+j])/2.0;
00668                   yy2[0]=(lookupt[loc3+j] -lookupt[loc5+j])/2.0; //derivative in y direction
00669                   yy2[1]=(lookupt[loc2+j] -lookupt[loc6+j])/2.0;
00670                   yy2[2]=(lookupt[loc9+j] -lookupt[loc1+j])/2.0;
00671                   yy2[3]=(lookupt[loc10+j]-lookupt[loc0+j])/2.0;
00672                   yy12[0]=(lookupt[loc2+j] -lookupt[loc6+j] -lookupt[loc11+j]+lookupt[loc12+j])/4.0; //cross-derivatives
00673                   yy12[1]=(lookupt[loc8+j] -lookupt[loc13+j]-lookupt[loc3+j] +lookupt[loc5+j] )/4.0;
00674                   yy12[2]=(lookupt[loc14+j]-lookupt[loc7+j] -lookupt[loc10+j]+lookupt[loc0+j] )/4.0;
00675                   yy12[3]=(lookupt[loc9+j] -lookupt[loc1+j] -lookupt[loc15+j]+lookupt[loc4+j] )/4.0;
00676                   bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00677                   u=yb;
00678                   t=xb;
00679                   poly[j]=0.0;
00680                   for (k = 3; k >= 0; k--) poly[j] = t*poly[j] + ((coeffs[k][3]*u+coeffs[k][2])*u + coeffs[k][1])*u + coeffs[k][0];
00681                   
00682                   //for 1st Fourier coefficient
00683                   if(poly[j] > vLCP && poly[j-1]   <= vLCP) indexL = j-1;
00684                 }
00685             }
00686               if(poly[i] > vRCP && poly[i-step]   <= vRCP)
00687             {
00688               for(j=i-step+1;j<=i;j++)
00689                 {
00690 
00691                   yy[0]=lookupt[loc0+j];
00692                   yy[1]=lookupt[loc1+j];
00693                   yy[2]=lookupt[loc2+j];
00694                   yy[3]=lookupt[loc3+j];
00695                   yy1[0]=(lookupt[loc1+j]-lookupt[loc4+j])/2.0; //derivative in x direction
00696                   yy1[1]=(lookupt[loc7+j]-lookupt[loc0+j])/2.0;
00697                   yy1[2]=(lookupt[loc8+j]-lookupt[loc3+j])/2.0;
00698                   yy1[3]=(lookupt[loc2+j]-lookupt[loc11+j])/2.0;
00699                   yy2[0]=(lookupt[loc3+j] -lookupt[loc5+j])/2.0; //derivative in y direction
00700                   yy2[1]=(lookupt[loc2+j] -lookupt[loc6+j])/2.0;
00701                   yy2[2]=(lookupt[loc9+j] -lookupt[loc1+j])/2.0;
00702                   yy2[3]=(lookupt[loc10+j]-lookupt[loc0+j])/2.0;
00703                   yy12[0]=(lookupt[loc2+j] -lookupt[loc6+j] -lookupt[loc11+j]+lookupt[loc12+j])/4.0; //cross-derivatives
00704                   yy12[1]=(lookupt[loc8+j] -lookupt[loc13+j]-lookupt[loc3+j] +lookupt[loc5+j] )/4.0;
00705                   yy12[2]=(lookupt[loc14+j]-lookupt[loc7+j] -lookupt[loc10+j]+lookupt[loc0+j] )/4.0;
00706                   yy12[3]=(lookupt[loc9+j] -lookupt[loc1+j] -lookupt[loc15+j]+lookupt[loc4+j] )/4.0;
00707                   bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00708                   u=yb;
00709                   t=xb;
00710                   poly[j]=0.0;
00711                   for (k = 3; k >= 0; k--) poly[j] = t*poly[j] + ((coeffs[k][3]*u+coeffs[k][2])*u + coeffs[k][1])*u + coeffs[k][0];
00712 
00713                   //for 1st Fourier coefficient
00714                   if(poly[j] > vRCP && poly[j-1]   <= vRCP) indexR = j-1;
00715                 }
00716             }
00717               if(poly2[i]> v2LCP && poly2[i-step] <= v2LCP)
00718             {
00719               for(j=i-step+1;j<=i;j++)
00720                 {
00721 
00722                   yy[0]=lookupt[loc0+j+ntest];
00723                   yy[1]=lookupt[loc1+j+ntest];
00724                   yy[2]=lookupt[loc2+j+ntest];
00725                   yy[3]=lookupt[loc3+j+ntest];
00726                   yy1[0]=(lookupt[loc1+j+ntest]-lookupt[loc4+j+ntest])/2.0; //derivative in x direction
00727                   yy1[1]=(lookupt[loc7+j+ntest]-lookupt[loc0+j+ntest])/2.0;
00728                   yy1[2]=(lookupt[loc8+j+ntest]-lookupt[loc3+j+ntest])/2.0;
00729                   yy1[3]=(lookupt[loc2+j+ntest]-lookupt[loc11+j+ntest])/2.0;
00730                   yy2[0]=(lookupt[loc3+j+ntest] -lookupt[loc5+j+ntest])/2.0; //derivative in y direction
00731                   yy2[1]=(lookupt[loc2+j+ntest] -lookupt[loc6+j+ntest])/2.0;
00732                   yy2[2]=(lookupt[loc9+j+ntest] -lookupt[loc1+j+ntest])/2.0;
00733                   yy2[3]=(lookupt[loc10+j+ntest]-lookupt[loc0+j+ntest])/2.0;
00734                   yy12[0]=(lookupt[loc2+j+ntest] -lookupt[loc6+j+ntest] -lookupt[loc11+j+ntest]+lookupt[loc12+j+ntest])/4.0; //cross-derivatives
00735                   yy12[1]=(lookupt[loc8+j+ntest] -lookupt[loc13+j+ntest]-lookupt[loc3+j+ntest] +lookupt[loc5+j+ntest] )/4.0;
00736                   yy12[2]=(lookupt[loc14+j+ntest]-lookupt[loc7+j+ntest] -lookupt[loc10+j+ntest]+lookupt[loc0+j+ntest] )/4.0;
00737                   yy12[3]=(lookupt[loc9+j+ntest] -lookupt[loc1+j+ntest] -lookupt[loc15+j+ntest]+lookupt[loc4+j+ntest] )/4.0;
00738                   bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00739                   u=yb;
00740                   t=xb;
00741                   poly2[j]=0.0;
00742                   for (k = 3; k >= 0; k--) poly2[j] = t*poly2[j] + ((coeffs[k][3]*u+coeffs[k][2])*u + coeffs[k][1])*u + coeffs[k][0];
00743 
00744                   //for 2nd Fourier coefficient
00745                   if(poly2[j] > v2LCP && poly2[j-1]   <= v2LCP) indexL2 = j-1;
00746                 }
00747             }
00748               if(poly2[i]> v2RCP && poly2[i-step] <= v2RCP)
00749             {
00750               for(j=i-step+1;j<=i;j++)
00751                 {
00752 
00753                   yy[0]=lookupt[loc0+j+ntest];
00754                   yy[1]=lookupt[loc1+j+ntest];
00755                   yy[2]=lookupt[loc2+j+ntest];
00756                   yy[3]=lookupt[loc3+j+ntest];
00757                   yy1[0]=(lookupt[loc1+j+ntest]-lookupt[loc4+j+ntest])/2.0; //derivative in x direction
00758                   yy1[1]=(lookupt[loc7+j+ntest]-lookupt[loc0+j+ntest])/2.0;
00759                   yy1[2]=(lookupt[loc8+j+ntest]-lookupt[loc3+j+ntest])/2.0;
00760                   yy1[3]=(lookupt[loc2+j+ntest]-lookupt[loc11+j+ntest])/2.0;
00761                   yy2[0]=(lookupt[loc3+j+ntest] -lookupt[loc5+j+ntest])/2.0; //derivative in y direction
00762                   yy2[1]=(lookupt[loc2+j+ntest] -lookupt[loc6+j+ntest])/2.0;
00763                   yy2[2]=(lookupt[loc9+j+ntest] -lookupt[loc1+j+ntest])/2.0;
00764                   yy2[3]=(lookupt[loc10+j+ntest]-lookupt[loc0+j+ntest])/2.0;
00765                   yy12[0]=(lookupt[loc2+j+ntest] -lookupt[loc6+j+ntest] -lookupt[loc11+j+ntest]+lookupt[loc12+j+ntest])/4.0; //cross-derivatives
00766                   yy12[1]=(lookupt[loc8+j+ntest] -lookupt[loc13+j+ntest]-lookupt[loc3+j+ntest] +lookupt[loc5+j+ntest] )/4.0;
00767                   yy12[2]=(lookupt[loc14+j+ntest]-lookupt[loc7+j+ntest] -lookupt[loc10+j+ntest]+lookupt[loc0+j+ntest] )/4.0;
00768                   yy12[3]=(lookupt[loc9+j+ntest] -lookupt[loc1+j+ntest] -lookupt[loc15+j+ntest]+lookupt[loc4+j+ntest] )/4.0;
00769                   bcucof(yy,yy1,yy2,yy12,1.0,1.0,coeffs);
00770                   u=yb;
00771                   t=xb;
00772                   poly2[j]=0.0;
00773                   for (k = 3; k >= 0; k--) poly2[j] = t*poly2[j] + ((coeffs[k][3]*u+coeffs[k][2])*u + coeffs[k][1])*u + coeffs[k][0];
00774 
00775                   //for 2nd Fourier coefficient
00776                   if(poly2[j] > v2RCP && poly2[j-1]   <= v2RCP) indexR2 = j-1;
00777                 } 
00778             }
00779             }
00780 
00781 
00782           //TO DEAL WITH SATURATION
00783           if(vLCP < minlookupt1)
00784             { 
00785               vLCP=vtest[0];
00786               indexL=ntest+1;
00787               //printf("take1 %f %f %f %f\n",vLCP,minlookupt1,vRCP,distance);
00788             }
00789           if(vLCP > maxlookupt1)
00790             {
00791               vLCP=vtest[ntest-1];
00792               indexL=ntest+1;
00793               //printf("take2 %f %f %f %f\n",vLCP,maxlookupt1,vRCP,distance);
00794             }
00795           if(vRCP < minlookupt1)
00796             { 
00797               vRCP=vtest[0];
00798               indexR=ntest+1;
00799               //printf("take3 %f %f %f %f\n",vRCP,minlookupt1,vLCP,distance);
00800             }
00801           if(vRCP > maxlookupt1)
00802             {
00803               vRCP=vtest[ntest-1];
00804               indexR=ntest+1;
00805               //printf("take4 %f %f %f %f\n",vRCP,maxlookupt1,vLCP,distance);
00806             }
00807           if(v2LCP < minlookupt2)
00808             { 
00809               v2LCP=vtest[0];
00810               indexL2=ntest+1;
00811             }
00812           if(v2LCP > maxlookupt2)
00813             {
00814               v2LCP=vtest[ntest-1];
00815               indexL2=ntest+1;
00816             }
00817           if(v2RCP < minlookupt2)
00818             { 
00819               v2RCP=vtest[0];
00820               indexR2=ntest+1;
00821             }
00822           if(v2RCP > maxlookupt2)
00823             {
00824               v2RCP=vtest[ntest-1];
00825               indexR2=ntest+1;
00826             }
00827 
00828           if(indexL == ntest+1 || indexR == ntest+1 || indexL2 == ntest+1 || indexR2 == ntest+1)
00829             {
00830               SATVALS2   += 1;
00831             }
00832           //else
00833           //{
00834 
00835           //We linearly interpolate in the look-up table for the 1st Fourier coefficient to retrieve the actual velocities
00836           if(indexL != ntest+1)  vLCP   = vtest[indexL] +(vLCP-(double)poly[indexL])   *(vtest[indexL+1] -vtest[indexL]) /((double)poly[indexL+1]  -(double)poly[indexL]  );
00837           if(indexR != ntest+1)  vRCP   = vtest[indexR] +(vRCP-(double)poly[indexR])   *(vtest[indexR+1] -vtest[indexR]) /((double)poly[indexR+1]  -(double)poly[indexR]  );
00838           //We linearly interpolate in the look-up table for the 2nd Fourier coefficient to retrieve the actual velocities
00839           if(indexL2 != ntest+1) v2LCP  = vtest[indexL2]+(v2LCP-(double)poly2[indexL2])*(vtest[indexL2+1]-vtest[indexL2])/((double)poly2[indexL2+1]-(double)poly2[indexL2]);
00840           if(indexR2 != ntest+1) v2RCP  = vtest[indexR2]+(v2RCP-(double)poly2[indexR2])*(vtest[indexR2+1]-vtest[indexR2])/((double)poly2[indexR2+1]-(double)poly2[indexR2]);
00841           
00842           /*-------------------------------------------------------------*/
00843           /* Calculation of observables                                  */
00844           /*-------------------------------------------------------------*/
00845           
00846           f1LCPc = f1LCPc*Kfourier;
00847           f1RCPc = f1RCPc*Kfourier;
00848           f1LCPs = f1LCPs*Kfourier;
00849           f1RCPs = f1RCPs*Kfourier;
00850           f2LCPc = f2LCPc*Kfourier;
00851           f2RCPc = f2RCPc*Kfourier;
00852           f2LCPs = f2LCPs*Kfourier;
00853           f2RCPs = f2RCPs*Kfourier;
00854           
00855           //We compute the Doppler velocity
00856           //lam0g[iii]  = (float)((vLCP+vRCP+v2LCP+v2RCP)/4.);//simple average. Need weights? REMINDER: SIGN CONVENTION: v<0 FOR MOTION TOWARD THE OBSERVER (BLUESHIFT)
00857           lam0g[iii]  = (float)((vLCP+vRCP)/2.);
00858           
00859           //We compute the l.o.s. magnetic field
00860           //B0g[iii]    = (float)((vLCP-vRCP+v2LCP-v2RCP)/2.0*magnetic);
00861           B0g[iii]    = (float)((vLCP-vRCP)*magnetic);
00862           
00863           //We compute the linewidth (in Angstroms)
00864           temp        = period/M_PI*sqrt(1.0/6.0*log((f1LCPc*f1LCPc+f1LCPs*f1LCPs)/(f2LCPc*f2LCPc+f2LCPs*f2LCPs)));
00865           tempbis     = period/M_PI*sqrt(1.0/6.0*log((f1RCPc*f1RCPc+f1RCPs*f1RCPs)/(f2RCPc*f2RCPc+f2RCPs*f2RCPs)));
00866           widthg[iii] = (float)((temp+tempbis)*sqrt(log(2.0)))*1000.; //we want the FWHM not the sigma of the Gaussian, in milliAngstoms          
00867           
00868           //We compute the linedepth
00869           //temp2       = period/2.0*sqrt(f1LCPc*f1LCPc+f1LCPs*f1LCPs)/sqrt(M_PI)/temp*exp(M_PI*M_PI*temp*temp/period/period);
00870           //temp2bis    = period/2.0*sqrt(f1RCPc*f1RCPc+f1RCPs*f1RCPs)/sqrt(M_PI)/tempbis*exp(M_PI*M_PI*tempbis*tempbis/period/period);
00871           //Idg[iii]    = (float)((temp2+temp2bis)/2.0);
00872           
00873           //We compute the linedepth (alternative formula)
00874           //temp2       = period/2.0/sqrt(M_PI)/temp   *pow((f1LCPc*f1LCPc+f1LCPs*f1LCPs),2./3.)/pow((f2LCPc*f2LCPc+f2LCPs*f2LCPs),1./6.);
00875           //temp2bis    = period/2.0/sqrt(M_PI)/tempbis*pow((f1RCPc*f1RCPc+f1RCPs*f1RCPs),2./3.)/pow((f2RCPc*f2RCPc+f2RCPs*f2RCPs),1./6.);
00876           temp2       = period/2.0*sqrt(f1LCPc*f1LCPc+f1LCPs*f1LCPs)/sqrt(M_PI)/FWHM*exp(M_PI*M_PI*FWHM*FWHM/period/period);  //to not use the second Fourier coefficient
00877           temp2bis    = period/2.0*sqrt(f1RCPc*f1RCPc+f1RCPs*f1RCPs)/sqrt(M_PI)/FWHM*exp(M_PI*M_PI*FWHM*FWHM/period/period);
00878           Idg[iii]    = (float)((temp2+temp2bis)/2.0);
00879           
00880           //We compute the continuum intensity
00881           //temp3       = (vLCP+v2LCP)/2.0/dv;
00882           //temp3bis    = (vRCP+v2RCP)/2.0/dv;
00883           //meanL=0.0;
00884           //meanR=0.0;
00885           //for(i=0;i<N;++i)
00886           //{
00887           //  meanL      += L[i];
00888           //      meanR      += R[i];
00889           //}
00890           //meanL=meanL/(double)N;
00891           //meanR=meanR/(double)N;
00892           //for(i=0;i<N;++i)
00893           //{
00894           //  meanL      += (temp2   /(double)N*exp(-(tune[i]-temp3)   *(tune[i]-temp3)   /temp   /temp));
00895           //  meanR      += (temp2bis/(double)N*exp(-(tune[i]-temp3bis)*(tune[i]-temp3bis)/tempbis/tempbis));
00896           //}
00897           //
00898           //I0g[iii]    = (float)((meanL+meanR)/2.0);
00899           
00900           temp3       = vLCP/dv;
00901           temp3bis    = vRCP/dv;
00902           meanL=0.0;
00903           meanR=0.0;
00904           for(i=0;i<N;++i)
00905             {
00906               meanL      += L[i];
00907               meanR      += R[i];
00908             }
00909           meanL=meanL/(double)N;
00910           meanR=meanR/(double)N;
00911           //minimumL=(L[0]+L[N-1])/2.*minimum/(double)N; //(L[0]+L[N-1])/2. estimate of the continuum
00912           //minimumR=(R[0]+R[N-1])/2.*minimum/(double)N;
00913           for(i=0;i<N;++i)
00914             {
00915               //meanL      += (minimumL*exp(-(tune[i]-temp3)   *(tune[i]-temp3)   /FWHM/FWHM));
00916               //meanR      += (minimumR*exp(-(tune[i]-temp3bis)*(tune[i]-temp3bis)/FWHM/FWHM));
00917               meanL      += (temp2   /(double)N*exp(-(tune[i]-temp3)   *(tune[i]-temp3)   /FWHM/FWHM));
00918               meanR      += (temp2bis/(double)N*exp(-(tune[i]-temp3bis)*(tune[i]-temp3bis)/FWHM/FWHM));
00919             }
00920           
00921           I0g[iii]    = (float)((meanL+meanR)/2.0);
00922           
00923           // }        
00924           
00925         }//if (distance <= Rsun)
00926       else
00927         {
00928           lam0g[iii]  = MISSINGRESULT;
00929           B0g[iii]    = MISSINGRESULT;
00930           widthg[iii] = MISSINGRESULT;
00931           Idg[iii]    = MISSINGRESULT;
00932           I0g[iii]    = MISSINGRESULT;
00933         }
00934       
00935     }//for iii
00936  }//end #pragma omp parallel
00937  
00938  
00939  free(cosi);
00940  free(sini);
00941  free(cos2i);
00942  free(sin2i);
00943  free(tune);
00944  free(angle);
00945  
00946  //fclose(fp);
00947  
00948  *MISSVALS=MISSVALS2;
00949  *SATVALS=SATVALS2;
00950  return error;
00951  
00952 }//end routine