proj/flatfield/apps/module_flatfield.c

/*
 * module_flatfield - Calculates rotational flatfields (first step) and produces cosmic ray records (first step)
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <jsoc_main.h>
#include <string.h>
#include <time.h>
#include <omp.h>
#include <module_flatfield.h>

char *module_name  = "module_flatfield";    //name of the module

#define kRecSetIn      "input_series" //name of the series containing the input filtergrams
#define kRecSetOut      "cosmic_ray_series" //name of the series containing the output series for cosmic rays // optional
#define kDSCosmic      "cosmic_rays" //name of the cosmic ray series
#define kDSFlatfield      "flatfield" 
#define cadence_name "cadence"  //cadence in sec
#define fid_name "fid" //name of fid argument
#define cam_name "camera" // 0 side, 1 front
#define fsnf_name "fsn_first"
#define fsnl_name "fsn_last"
#define datumn  "datum"


#define minval(x,y) (((x) < (y)) ? (x) : (y))                                        
#define maxval(x,y) (((x) < (y)) ? (y) : (x))

ModuleArgs_t module_args[] = {
     {ARG_STRING, kRecSetIn, "",  "Input data series."},
     {ARG_STRING, datumn, "", "datum string"},
     {ARG_INT, kDSCosmic, "", "Cosmic rays flag"},
     {ARG_INT, kDSFlatfield, "", "Flatfield flag"},
     {ARG_INT, cadence_name, "", "Cadence in sec"},
     {ARG_INT, fid_name, "", "FID"},
     {ARG_INT, cam_name, "", "Camera"},
     {ARG_INT, fsnf_name, "0"},
     {ARG_INT, fsnl_name, "2147483647"},
  {ARG_FLOAT, "rmax_crfix", "0.98",
      "Maximum value of r/R for cosmic ray detection"},
     {ARG_STRING, kRecSetOut, "default",  "Cosmic ray output series."},
     
     {ARG_END}
};









long sign(double);
double mat_rh(long[], double[], double[], int);
void nine_diag(long[], double[], int, int, double[], long[], double);
void blockiter(double[], double[], double[], double*, double[], int, int, double);
void tridag(double veca[], double vecb[], double vecc[], double vecr[], double vecu[], int n);
void mat_mult(double[], double[], double[], double[], double[], int);
void printtime();
void highpass(int, int, double, double[nx][ny]);
void derotation(double, double, double, double, double, double, double [2][nx][ny]);

void limb_darkening(double radius, double cent_x, double cent_y, double *b, int order, double *limb_dark);
int flatfield(double *rhsp, double *rhsm,short *badpix, int pairs, double *gout, double *param, struct code_param, double deltat);
void apod_circ(float rad, float nb, float offx, float offy, float *vd);


int get_flat(char *query, int camera, float *flatfield, float *offpoint, 
         short *bad);

int get_latest_bad(TIME t_0, int camera, const char *fname, const char *segname, short *bad, long long *recnum);

int get_latest_flat(TIME t_0, int camera, const char *fname, const char *segname, float *offpoint, long long *recnum, int *focus);



int write_flatfield_fid(DRMS_Array_t *arrout_new, int camera, long long recnum_offpoint, 
                TIME t_0, int focus,int fid,
            struct rotpar rot_new);
/*-------------------------------------------------------------*/
/*                                                             */
/*   DoIt is the entry point of the module                     */
/*   the name MUST be DoIt for a DRMS module                   */
/*                                                             */
/*-------------------------------------------------------------*/

int DoIt(void)
{

#include "module_flatfield_const.h"

  //read out input parameters


  const char *inRecQuery = cmdparams_get_str(&cmdparams, kRecSetIn, NULL); //cmdparams is defined in jsoc_main.h
  const char *cosmic_ray_series = cmdparams_get_str(&cmdparams, kRecSetOut, NULL); //cmdparams is defined in jsoc_main.h
  if (strcmp(cosmic_ray_series, "default") == 0){cosmic_ray_series=filename_cosmic; }

  int cosmic_flag=cmdparams_get_int(&cmdparams, kDSCosmic, NULL);
  int flatfield_flag=cmdparams_get_int(&cmdparams, kDSFlatfield, NULL);
  int cadence=cmdparams_get_int(&cmdparams, cadence_name, NULL);

  int fid=cmdparams_get_int(&cmdparams, fid_name, NULL);
  int cam=cmdparams_get_int(&cmdparams, cam_name, NULL)-1;
  int fsn_first=cmdparams_get_int(&cmdparams, fsnf_name, NULL);
  int fsn_last=cmdparams_get_int(&cmdparams, fsnl_name, NULL);
  float crrmax = cmdparams_get_float (&cmdparams, "rmax_crfix", NULL);

  const char *datum=cmdparams_get_str(&cmdparams, datumn, NULL);
 

  if (fid < minfid || fid > maxfid){printf("Not an observable FID\n"); return 1;}

  int  status, status1, status2; 
  int status0=DRMS_SUCCESS;
  int stat = DRMS_SUCCESS; 
  int status_flatfield=DRMS_SUCCESS;
  int status_flatfield_relative=DRMS_SUCCESS;
 

  DRMS_RecordSet_t *data, *data_flat, *data_last, *dataout;
  DRMS_Record_t  *ff_last=NULL;
  DRMS_Array_t *arrin0;
  double *arrinL0; //pointer to data
  DRMS_Segment_t *segin = NULL, *segin_flat=NULL, *segout = NULL, *segin_offpoint=NULL, *segout_val=NULL, *segout_sig=NULL;
  DRMS_Record_t *rec_flat, *recout, *reclink_off;

  DRMS_Type_t type_time = DRMS_TYPE_TIME;
  DRMS_Type_t type_int = DRMS_TYPE_INT;
  DRMS_Type_t type_float = DRMS_TYPE_FLOAT;
  DRMS_Type_t type_double = DRMS_TYPE_DOUBLE;

  // int axisin[2]={nx,ny};                   
                                 //size of input arrays

  double deltat=(double)cadence;
  printf("cadence %lf seconds\n", deltat);

  TIME t_0, t_stamp;
  int focus,camid,fsns;
  
  TIME   interntime=sscan_time("2010.03.18_22:12:17.77_UTC");
  int axisbad[1];


  int i,j,k,km1,c,ki,iii,jjj,kkk;                                                   //loop variables
  int signid=0;
  int nfr;
 
  /********************************************************************************************************************************/
  //define and initialize quantities
  /********************************************************************************************************************************/



  //sigma for each FID (cosmic rays)

 
  //  double *sigma=(double *)(malloc(nx*ny*sizeof(double)));
  float *limw=(float *)(malloc(nx*ny*sizeof(float)));


  //DRMS_arrays for output arrays
  DRMS_Array_t *arr_flat, *arr_offpoint;
  DRMS_Array_t *arrout_new;
  DRMS_Array_t *arrout, *arrout_val, *arrout_sig;

  int *cosmic_ray_data;
  float *val_data, *sig_data;

  //Array pointers for storage
 
  short *bad;
  short *badpix;
  double *flati;
 
  float *apod;

  float *flat_off, *flat_relative;
  float *app_flat;

  app_flat=(float *)(malloc(nx*ny*sizeof(float)));
  flat_off=(float *)(malloc(nx*ny*sizeof(float)));
  flat_relative=(float *)(malloc(nx*ny*sizeof(float)));
  badpix=(short *)(malloc(nx*ny*sizeof(short)));
  flati=(double *)(malloc(nx*ny*sizeof(double)));
 
  apod=(float *)(malloc(nx*ny*sizeof(float)));

  float *flatfield_new;
  float *fflat;
  
 

 

  double flat[nx][ny];
  for (j=0; j<ny; ++j) for (i=0; i<nx; ++i){flat[i][j]=0.0;}
 


  double *rhsp, *rhsm;
  double *limb_dark;
  int *count_p, *count_m;

  int count, ccount,ccount1,ccount2, count_filtergram; // counters 
  int count_flatfields;

  int last, current, next, holder;
  int index_last, index_current, index_next;
  float wlr;
  float *rr, *a1,*a2, *w1,*w2,*cmc;
  int flatfound, flatrec_index;

  rr=(float *)(malloc(nx*ny*sizeof(float)));
  a1=(float *)(malloc(nx*ny*sizeof(float)));
  a2=(float *)(malloc(nx*ny*sizeof(float)));
  w1=(float *)(malloc(nx*ny*sizeof(float)));
  w2=(float *)(malloc(nx*ny*sizeof(float)));
  cmc=(float *)(malloc(nx*ny*sizeof(float)));
  

  rhsp=(double *)(malloc(nx*ny*sizeof(double)));
  rhsm=(double *)(malloc(nx*ny*sizeof(double)));
         
  count_p=(int *)(malloc(nx*ny*sizeof(int)));
  count_m=(int *)(malloc(nx*ny*sizeof(int)));

  limb_dark=(double *)(malloc(nx*ny*sizeof(double)));

  //cosmic ray
  
  
  float arrimg[3][ny][nx];
  int *cosmicarr;
  cosmicarr=(int *)(malloc(nx*ny*sizeof(int)));
  float *cosmicval;
  cosmicval=(float *)(malloc(nx*ny*sizeof(float)));
  float *cosmicsig;
  cosmicsig=(float *)(malloc(nx*ny*sizeof(float)));
  float diff_forward, diff_backward;

  int fsnarr[3]={0 , 0, 0};

  unsigned char *cmarr=(unsigned char *)(malloc(nx*ny*sizeof(unsigned char)));

  char *query_flat, *query_flat_ref, *query_flat_relative;
  char flat_default[256]={"none"};
 
  double *param;
  param=(double *)(malloc(6*sizeof(double)));

  
  //parameters to identify records
  TIME tobs_link[2];
  long long recnum_offpoint, recnum_badpix;
  struct rotpar rot_cur;
  struct rotpar rot_new;

  int rotpairs;


  printf("START!\n");

  /***********************************************************************************************************/
  /*CHECK WHETHER THE FLATFIELD OUTPUT SERIES EXIST                                                                    */
  /***********************************************************************************************************/
    
       drms_series_exists(drms_env, filename_flatfield_fid, &status);
      if (status == DRMS_ERROR_UNKNOWNSERIES)
    {
      printf("Output series %s doesn't exist\n",filename_flatfield_fid);       //if the output series does not exit
      return 1;                                        //we exit the program
    } 
      if (status == DRMS_SUCCESS)
    {
      printf("Output series %s exists.\n",filename_flatfield_fid);
    }

  /***********************************************************************************************************/
  /*CHECK WHETHER THE COSMIC RAY OUTPUT SERIES EXIST                                                                    */
  /***********************************************************************************************************/
    
      drms_series_exists(drms_env, cosmic_ray_series, &status);
      if (status == DRMS_ERROR_UNKNOWNSERIES)
    {
      printf("Cosmic ray series %s doesn't exist\n",cosmic_ray_series);       //if the output series does not exit
      return 1;                                        //we exit the program
    } 
      if (status == DRMS_SUCCESS)
    {
      printf("Cosmic ray series %s exists.\n",cosmic_ray_series);
    }
  /**************************************************************************************************************/


      /**********************************/
      //parallelization
      /**********************************/


      int nthreads; 
      nthreads=omp_get_max_threads();
      //nthreads=omp_get_num_procs();                                      //number of threads supported by the machine where the code is running
      //omp_set_num_threads(nthreads);                                     //set the number of threads to the maximum value
      printf("Number of threads run in parallel = %d \n",nthreads);

  
  /**************************************************************************************************************/
  //built query strings
  /**************************************************************************************************************/
            
 char timefirst[256]="";
 strcat(timefirst, datum);
 strcat(timefirst, "_00:00:00.00_TAI");

 char timelast[256]="";
 strcat(timelast, datum);
 strcat(timelast, "_23:59:59.99_TAI");

 int pad;
 pad=(int)(time_limit/filtergram_cadence+1.0);


 TIME tfirst, tlast;

 char query0[256]="";

if (fsn_first == 0 || fsn_last == 2147483647)
  {

    tfirst=sscan_time(timefirst);
    tlast=sscan_time(timelast);

    char datefirst[256]="";
    sprint_ut(datefirst, tfirst-time_limit);

    char datelast[256]="";
    sprint_ut(datelast, tlast+time_limit);


 strcat(query0, inRecQuery);
 strcat(query0, "[");
 strcat(query0, datefirst);
 strcat(query0, "-");
 strcat(query0, datelast);
 strcat(query0, "][?FID=");
 char ffnumb[2]={""};
 sprintf(ffnumb, "%5.5d", fid);
 strcat(query0, ffnumb);
 strcat(query0, "?][?CAMERA=");
 char fnumb[2]={""};
 sprintf(fnumb, "%1.1d", cam+1);
 strcat(query0, fnumb);
 strcat(query0, "?]");

  printf("query string: %s\n", query0);

  }



  if (fsn_first != 0 && fsn_last != 2147483647)
    {

      tfirst=0;
      tlast=3881520000.0;

      strcat(query0, inRecQuery);
      strcat(query0, "[][");
       char fsnf[10]={""};
       sprintf(fsnf, "%d", fsn_first-pad);
       strcat(query0, fsnf);
      strcat(query0, "-");
      char fsnl[10]={""};
      sprintf(fsnl, "%d", fsn_last+pad);
       strcat(query0, fsnl);
      strcat(query0, "][?FID=");
      char ffnumb[2]={""};
     sprintf(ffnumb, "%5.5d", fid);
      strcat(query0, ffnumb);
      strcat(query0, "?][?CAMERA=");
      char fnumb[2]={""};
     sprintf(fnumb, "%1.1d", cam+1);
     strcat(query0, fnumb);
      strcat(query0, "?]");
      printf("query string: %s\n", query0);  
    }

  /************************************************************************/


  
  /*************************************************************************/
  //open records and stage data
  /*************************************************************************/


  
      data     = drms_open_records(drms_env,query0,&stat);
      if (data == NULL || stat != 0 || data->n == 0){printf("can not open records\n"); return 1;}

      drms_stage_records(data, 1, 0); // stage needed records from tape and wait until ready.


      /************************************************************************/
      //define keyword arrays
      /************************************************************************/


  int nRecs=data->n;
  printf("number of records %d\n", nRecs);




  DRMS_Record_t *rec0[nRecs];
 
  //array for keyvalues
  int keyvalue_cam[nRecs];
  int keyvalue_wl[nRecs];
  int keyvalue_pl[nRecs];
  DRMS_Type_Value_t  keyvalue_time;
  long tmind[nRecs+1];
  int keyvalue_fid[nRecs];
  char *keyvalue_iss[nRecs];
  char *keyvalue_flatnumb[nRecs];
  int keyvalue_flatr[nRecs];
  int keyvalue_fsn[nRecs];
  double time_fl[nRecs];
  //int cosmic_ray_check[nRecs];

  //arrays for Lev 1 solar parameters   
  float keyvalue_rsun[nRecs];
  double keyvalue_dsun_obs[nRecs];
  float keyvalue_p0[nRecs];
  float keyvalue_b0[nRecs];
  float keyvalue_X0[nRecs];
  float keyvalue_Y0[nRecs];
  float keyvalue_vrad[nRecs];

  

  int present[nRecs];        //indices
  int present_forward[nRecs];
  int present_backward[nRecs];
  int index[nRecs];
  int statarr[nRecs];

  /********************************************************************************/


  /****************************************************/
  //get time reference
  /****************************************************/

   
 
  t_0 = drms_getkey_time(data->records[nRecs/2],keytobs,&status);
  fsns=drms_getkey_int(data->records[nRecs/2],keyfsn,&status);
  focus = drms_getkey_int(data->records[nRecs/2],keyfocus,&status);
  
  char tmstr[256]="";
  strcat(tmstr, datum);
  strcat(tmstr, "_00:00:00.00_TAI");

  t_stamp=sscan_time(tmstr)+24.0*60.0*60.0; 

 
  /************************************************************/


  int update_stat=2;   // 2: not enough data // 1: no update required // 0: flatfield updated 
    

  /*************************************************************/
  //do camera test
  /************************************************************/

  if (stat == DRMS_SUCCESS && data != NULL && nRecs > 0)
    {

      nfr=0;
      for (k=0; k<nRecs; ++k)
    {
      statarr[k]=0;
      rec0[k]=data->records[k];
      keyvalue_cam[k]=drms_getkey_int(rec0[k],keycam,&status);
      if (status != 0 || (keyvalue_cam[k] != cam_id_front && keyvalue_cam[k] != cam_id_side)) statarr[k]=statarr[k]+1024;
      ++nfr;
    }

     
     
   
     
      
   

  
      

      /********************************************************/
      //read records and test keyword range
      /********************************************************/
     

      for (k=0; k<nRecs; ++k){
        keyvalue_time=drms_getkey(rec0[k], keytobs, &type_time, &status);       //read "T_OBS" of filtergram
        time_fl[k] = keyvalue_time.time_val;
        tmind[k]=(long)time_fl[k]-(long)interntime;
    
      
                     //tmind: time since first filtergram (in s) 

    keyvalue_fsn[k]=drms_getkey_int(rec0[k], keyfsn, &status);

    keyvalue_wl[k]=drms_getkey_int(rec0[k],keywl,&status);               //wavelength ID
    keyvalue_pl[k]=drms_getkey_int(rec0[k],keypl,&status);       // Polarization ID
      

    //test level_1 keywords


    keyvalue_rsun[k]=drms_getkey_float(rec0[k], lev1_r_sun,&status); // /drms_getkey_float(rec0[k], lev1_imsc,&status0);
    if (keyvalue_rsun[k] < rsun_min || keyvalue_rsun[k] > rsun_max ||  isnan(keyvalue_rsun[k]))statarr[k]=statarr[k]+32; 

    keyvalue_dsun_obs[k]=drms_getkey_double(rec0[k], lev1_dist ,&status);
    if (status !=0 || keyvalue_dsun_obs[k] < dsun_obs_min || keyvalue_dsun_obs[k] > dsun_obs_max || isnan(keyvalue_dsun_obs[k])) statarr[k]=statarr[k]+16;

    keyvalue_p0[k]=drms_getkey_float(rec0[k], lev1_p0, &status); 
    if (status !=0 || keyvalue_p0[k] < p0_min || keyvalue_p0[k] > p0_max || isnan(keyvalue_p0[k])) statarr[k]=statarr[k]+4;

    keyvalue_b0[k]=drms_getkey_float(rec0[k], lev1_b0, &status);
    if (status !=0 || keyvalue_b0[k] < b0_min || keyvalue_b0[k] > b0_max ||  isnan(keyvalue_b0[k])) statarr[k]=statarr[k]+8;

    keyvalue_X0[k]=drms_getkey_float(rec0[k], lev1_x0, &status1);
    keyvalue_Y0[k]=drms_getkey_float(rec0[k], lev1_y0,&status2);
    if (status1 !=0  || status2 != 0 || keyvalue_X0[k] < X0_min || keyvalue_X0[k] > X0_max || isnan(keyvalue_X0[k]) || keyvalue_Y0[k] < Y0_min || keyvalue_Y0[k] > Y0_max || isnan(keyvalue_Y0[k])) statarr[k]=statarr[k]+64;


    keyvalue_fid[k]=drms_getkey_int(rec0[k],fidkey,&status);
    if (status !=0 || keyvalue_fid[k] < minfid || keyvalue_fid[k] > maxfid || isnan(keyvalue_fid[k])) statarr[k]=statarr[k]+256;
    keyvalue_iss[k]=drms_getkey_string(rec0[k],isskey,&status);
    if (status !=0 || strcmp(keyvalue_iss[k], "CLOSED") != 0) statarr[k]=statarr[k]+1;

    keyvalue_flatnumb[k]=drms_getkey_string(rec0[k],flatnkey,&status); 

    
    keyvalue_vrad[k]=drms_getkey_float(rec0[k], lev1_vr, &status);
    if (keyvalue_vrad[k] < vrad_min || keyvalue_vrad[k] > vrad_max) statarr[k]=statarr[k]+512;

   

      }
  

      /******************************/ 
      //calculate average for solar parameters /
      /******************************/


      float R_SUN=0.0;
      float XX0=0.0;
      float YY0=0.0;
      float P_ANG=0.0;
      float B_ANG=0.0;
      float dist=0.0;

      float dcount=0.0;
      for (k=0; k<nRecs; ++k){
    //printf("rsun %f %f %f %f %f %lf %d %d\n", keyvalue_rsun[k], keyvalue_X0[k], keyvalue_Y0[k], keyvalue_p0[k], keyvalue_b0[k], keyvalue_dsun_obs[k], statarr[k], keyvalue_fid[k]);
    if (statarr[k] == 0){

    R_SUN=R_SUN+keyvalue_rsun[k];
    XX0=XX0+keyvalue_X0[k];
    YY0=YY0+keyvalue_Y0[k];
    P_ANG=P_ANG-keyvalue_p0[k]/180.0f*M_PI;
    B_ANG=B_ANG+keyvalue_b0[k]/180.0f*M_PI;
    dist=dist+(float)keyvalue_dsun_obs[k]/oneau;
    dcount=dcount+1.0f;
    }
      }
    


      //
  
      if (dcount > 0.0){
    R_SUN=R_SUN/dcount;
    XX0=XX0/dcount;
    YY0=YY0/dcount;
    P_ANG=P_ANG/dcount;
    B_ANG=B_ANG/dcount;
    dist=dist/dcount;
     
       *(param+0)=(double)R_SUN;  //R_SUN in pixel
       *(param+1)=(double)XX0;    //center x in pixel
       *(param+2)=(double)YY0;    //center y in pixel
       *(param+3)=(double)P_ANG;  //P-angle in rad // 
       *(param+4)=(double)B_ANG;  //B-angle in rad
       *(param+5)=(double)dist;   //distance in AU
      } 
      else
    {
      printf("Invalid lev1 keywords\n");
      return 1;
    }
    
     
      printf("disk center and radius: %f %f %f\n", XX0, YY0,R_SUN);         ;
  

      /**************************************************************/
      //Start flatfield calculations
      /**************************************************************/

      /************************/
      //get valid frames for flatfield calculation
      /************************/

      if (cam ==0)camid=cam_id_side;
      if (cam ==1)camid=cam_id_front;



      count_flatfields=0;

      for (j=0; j<ny; ++j) for (i=0; i<nx; ++i){rhsm[j*nx+i]=0.0; rhsp[j*nx+i]=0.0; count_p[j*nx+i]=0; count_m[j*nx+i]=0; flati[j*nx+i]=0.0; cmarr[j*nx+i]=1;}

    
    for (i=0; i<nRecs; ++i){present_backward[i]=0; present_forward[i]=0; present[i]=0;}
    //for (i=0; i<nRecs; ++i){present[i]=0; index[i]=-1;}

      //count filtergrams, and number of valid pairs
      //count_filtergram=0; 
   
      for (k=0; k<nRecs; ++k){present[k]=1; index[k]=k;}  //loop over all filtergrams: search for valid filtergrams
      printf("number of filtergrams %d\n", nRecs);
      dataout = drms_create_records(drms_env,nRecs,(char *)cosmic_ray_series,DRMS_PERMANENT,&stat); //create ALL cosmic_ray records

    //check each forward and backward pair of images for identity
    //check for wavelength, polarization,  P-angle, center, and solar radius identity, and for time difference being nominal cadence

    
      count_filtergram=nRecs;

      //if (count_filtergram >=3){  //if at least 3 filtergrams of this type

      
      printf("filtergram with FID  %d: %d\n", fid, count_filtergram);
    for (c=1; c<count_filtergram; ++c){

//backward pairs
      present_backward[c]=present[c];
      if (statarr[c] != 0 || statarr[c-1] != 0 || keyvalue_wl[c] != keyvalue_wl[c-1] || keyvalue_pl[c] != keyvalue_pl[c-1] || (tmind[c]-tmind[c-1]) != (long)deltat || fabs(keyvalue_p0[c]- keyvalue_p0[c-1]) > 0.2 || sqrt(pow(keyvalue_X0[c]-keyvalue_X0[c-1],2)+pow(keyvalue_Y0[c]-keyvalue_Y0[c-1],2)) > limit_centerdiff[cam] || fabs(keyvalue_rsun[c]-keyvalue_rsun[c-1]) > limit_rsundiff){present_backward[c]=0;}

      if (sqrt(pow(keyvalue_X0[c]-keyvalue_X0[c-1],2)+pow(keyvalue_Y0[c]-keyvalue_Y0[c-1],2)) > limit_centerdiff_cosmic) statarr[c]=statarr[c]+128;
      //  printf("%d %d %ld %f %f\n", keyvalue_wl[index[c]]-keyvalue_wl[index[c-1]],keyvalue_pl[index[c]]- keyvalue_pl[index[c-1]], (tmind[index[c]]-tmind[index[c-1]]), fabs(keyvalue_p0[index[c]]- keyvalue_p0[index[c-1]]), fabs(keyvalue_rsun[index[c]]-keyvalue_rsun[index[c-1]]));
    }


    //forward pairs
    ccount=count_filtergram-1; //number of pairs (forward or backward)
    for (c=0; c<(count_filtergram-1); ++c){
      present_forward[c]=present[c];
      if (statarr[c] != 0 || statarr[c+1] != 0 || keyvalue_wl[c] != keyvalue_wl[c+1] || keyvalue_pl[c] != keyvalue_pl[c+1] || (tmind[c+1]-tmind[c]) != (long)deltat || fabs(keyvalue_p0[c+1]- keyvalue_p0[c]) > 0.2 || sqrt(pow(keyvalue_X0[c]-keyvalue_X0[c+1],2)+pow(keyvalue_Y0[c]-keyvalue_Y0[c+1],2)) > limit_centerdiff[cam] || fabs(keyvalue_rsun[c]-keyvalue_rsun[c+1]) > limit_rsundiff){present_forward[c]=0; --ccount;}

      if (sqrt(pow(keyvalue_X0[c]-keyvalue_X0[c+1],2)+pow(keyvalue_Y0[c]-keyvalue_Y0[c+1],2)) > limit_centerdiff_cosmic || fabs(keyvalue_X0[c]-XX0) > limit_offpoint || fabs(keyvalue_Y0[c]-YY0) > limit_offpoint) statarr[c]=statarr[c]+128;
    }

    printf("number of filtergrams, valid pairs of frames %d \t %d \n", count_filtergram, ccount);



    //  for (k=0; k < count_filtergram; ++k) printf("%d \t %d \t %s \t %d\t %d \t %d \t %ld   \t %d  \t %f \t %f \t  %f %f %d %d %d %d\n", k, keyvalue_fsn[k], keyvalue_iss[k], keyvalue_wl[k], keyvalue_pl[k], keyvalue_cam[k], tmind[k], keyvalue_fid[k], keyvalue_rsun[k], keyvalue_X0[k], keyvalue_Y0[k], keyvalue_vrad[k], statarr[k], present[k], present_forward[k], present_backward[k]);



    /*******************************************/
    /********************************************/

    status_flatfield_relative=get_latest_flat(t_stamp, cam+1, filename_offpoint, segmentname_offpoint, flat_relative, &recnum_offpoint, &focus);




     if (status_flatfield_relative == 0)
       {
         printf("Relative flatfield read: record number:%ld\n", recnum_offpoint);
       }
     else 
       {
         printf("WARNING: Could not read relative flatfield");
       }



     
      


     /******************************************/
     //calculate limit parameters for cosmic ray detection
     /******************************************/


#pragma omp parallel for private(jjj,iii)      
          for (jjj=0; jjj<ny; ++jjj) for (iii=0; iii<nx; ++iii)
        {
                  
          rr[jjj*nx+iii]=sqrt(((float)iii-XX0)*((float)iii-XX0)+((float)jjj-YY0)*((float)jjj-YY0))/R_SUN;
          a1[jjj*nx+iii]=coef0[cam][0]+coef0[cam][1]*rr[jjj*nx+iii]+coef0[cam][2]*rr[jjj*nx+iii]*rr[jjj*nx+iii];
          a2[jjj*nx+iii]=coef1[cam][0]+coef1[cam][1]*rr[jjj*nx+iii]+coef1[cam][2]*rr[jjj*nx+iii]*rr[jjj*nx+iii];
          w1[jjj*nx+iii]=coef2[cam][0]+coef2[cam][1]*rr[jjj*nx+iii]+coef2[cam][2]*rr[jjj*nx+iii]*rr[jjj*nx+iii];
          w2[jjj*nx+iii]=w1[jjj*nx+iii];
          cmc[jjj*nx+iii]=coef4[cam][0]+coef4[cam][1]*rr[jjj*nx+iii]+coef4[cam][2]*rr[jjj*nx+iii]*rr[jjj*nx+iii];

        }







      /**********************************************************************************/
      //loop over filtergrams: obtain cosmic ray records and get input arrays for rotational flatfield calculations
      /**********************************************************************************/

      printf("cosmic ray detection\n");

      ccount1=0;
      ccount2=0;
      last=2; current=0; next=1;
      double avg=0.0;
      query_flat_ref=flat_default;
      printtime();

      printf("count filtergram %d\n", count_filtergram);

      for (kkk=0; kkk<(count_filtergram+1); ++kkk){ // loop starts at 0
      
        if (cosmic_flag || ccount >= cthreshold[cam])
          {

        km1=kkk-1;

        /**************************/
        //read in filtergram

        if (kkk < count_filtergram)
          {
            printf("fsn %d %d\n", kkk, keyvalue_fsn[kkk]);
     
            segin    = drms_segment_lookupnum(rec0[kkk], 0);
            arrin0= drms_segment_read(segin, type_double, &status);
        

      
            if (status != DRMS_SUCCESS || arrin0 == NULL)
              {
        
            printf("Error: there is a problem with the filtergram number %d \n", kkk); //if there is a problem with the filtergram
            --ccount;

            present_backward[kkk]=0;
            if (kkk < (count_filtergram-1)) present_backward[kkk+1]=0;

            present_forward[kkk]=0;
            if (kkk > 0 ) present_forward[kkk-1]=0;

            present[kkk]=0;
            for (jjj=0; jjj<ny; ++jjj) for (iii=0; iii<nx; ++iii)   arrimg[next][jjj][iii]=NAN;
            flatfound=1;
              }
            else

              {
         
            printf("filtergram with FSN %d read\n",keyvalue_fsn[kkk]);
            arrinL0 = arrin0->data;


            /**************************************/
            //flatfield query

            query_flat =keyvalue_flatnumb[kkk];
            if (strcmp(query_flat, query_flat_ref) != 0)
              {
                printf("%d %d read new flatfield\n", k, keyvalue_fsn[kkk]);

         
                status_flatfield=get_flat(query_flat, cam+1, app_flat, flat_off,  badpix);
                flatfound=status_flatfield;

                if (status_flatfield ==0)
                  {
                printf("Flatfield %s read\n", query_flat);

                  }
                else

                  {
                printf("Failure reading flatfield / Use unity flatfield\n");

                status=get_latest_bad(t_stamp, cam+1, filename_badpix, segmentname_badpix, badpix, &recnum_badpix); //retrieve bad for links a
                if (status != 0){printf("could not find bad pixel list"); return 1;}
                else {printf("recnum badpix %ld\n", recnum_badpix);}
                  }
         


                query_flat_ref=query_flat;
              }

            if (status_flatfield != 0)
              {
                --ccount;
                statarr[kkk]=statarr[kkk]+2;

                present_forward[kkk]=0;
                if (kkk > 0 ) present_forward[kkk-1]=0;

                present_backward[kkk]=0;
                if (kkk < (count_filtergram-1)) present_backward[kkk+1]=0;
                printf("no flatfield for filtergram %d\n", kkk);
              }
    
         
            if (flatfound == 0) for (jjj=0; jjj<ny; ++jjj) for (iii=0; iii<nx; ++iii)   arrimg[next][jjj][iii]=(float)arrinL0[jjj*nx+iii]*app_flat[jjj*nx+iii]/flat_relative[jjj*nx+iii];           
            if (flatfound != 0) for (jjj=0; jjj<ny; ++jjj) for (iii=0; iii<nx; ++iii)   arrimg[next][jjj][iii]=(float)arrinL0[jjj*nx+iii];

            fsnarr[next]=keyvalue_fsn[kkk];

            keyvalue_flatr[kkk]=flatfound;
              }
            drms_free_array(arrin0);
          }


        
        count=-1;

        /**********************************************/
        //cosmic ray detection

        if (kkk >= 2 && kkk < count_filtergram)
          {
            index_last=kkk-2;
            index_current=kkk-1;
            index_next=kkk;

            if (statarr[index_current] < 4)
              {     

            flatrec_index=keyvalue_flatr[index_last]+keyvalue_flatr[index_current]+keyvalue_flatr[index_next]; //==0 if all three flatfields could be retrieved, ==3 if none has been retrieved

    
    
            if (statarr[index_last] < 4 && statarr[index_next] < 4 && (tmind[index_current]-tmind[index_last]) < time_limit && (tmind[index_next]-tmind[index_current]) < time_limit && (flatrec_index == 0 || flatrec_index == 3))
              {
                count=0;
#pragma omp parallel 
                {
#pragma omp for private(jjj,iii,wlr)       
                  for (jjj=0; jjj<ny; ++jjj) for (iii=0; iii<nx; ++iii)
                {
    
//                if (rr[jjj*nx+iii] < rad_cosmic_ray)
                  if (rr[jjj*nx+iii] < crrmax)
                    {
                      wlr=((float)(keyvalue_wl[km1] % 20)-5.0)/2.0*lambda_sep-keyvalue_vrad[km1]/v_c*lambda0-(cos(keyvalue_p0[km1]/180.*M_PI)*((float)iii-keyvalue_X0[km1])-sin(keyvalue_p0[km1]/180.*M_PI)*((float)jjj-keyvalue_Y0[km1]))*cos(keyvalue_b0[km1]/180.0*M_PI)/keyvalue_rsun[km1]*cpa.rotcoef0*radsun_mm/v_c*lambda0;
                      limw[jjj*nx+iii]=a1[jjj*nx+iii]*exp(-(wlr-coef3[cam][0])*(wlr-coef3[cam][0])/2.0/w1[jjj*nx+iii]/w1[jjj*nx+iii])+a2[jjj*nx+iii]*exp(-(wlr-coef3[cam][1])*(wlr-coef3[cam][1])/2.0/w2[jjj*nx+iii]/w2[jjj*nx+iii])+cmc[jjj*nx+iii];
                    }
                }


    
            
      //head= NULL;
#pragma omp for private(iii,jjj,diff_forward,diff_backward)
                  for (jjj=0; jjj<ny; ++jjj) for (iii=0; iii<nx; ++iii)
                {
                  cmarr[jjj*nx+iii]=1;
       
//                if (rr[jjj*nx+iii] < rad_cosmic_ray)
                  if (rr[jjj*nx+iii] < crrmax)
                    {
                      diff_backward=arrimg[current][jjj][iii]-arrimg[last][jjj][iii];
                      diff_forward=arrimg[current][jjj][iii]-arrimg[next][jjj][iii];

                      if (diff_forward < 0.0 && diff_backward < 0.0) signid=1; else signid=0;

                      if (fabs(diff_forward) > (factor[cam][signid]*limw[jjj*nx+iii]) && fabs(diff_backward) > (factor[cam][signid]*limw[jjj*nx+iii])){
                    
                    if (badpix[jjj*nx+iii]){
                      if (!isnan(arrimg[last][jjj][iii]) && !isnan(arrimg[current][jjj][iii]) && !isnan(arrimg[next][jjj][iii])){
#pragma omp critical(detect)
                        {
                          cosmicarr[count]=jjj*nx+iii;
                          if (diff_forward > 0.0 && diff_backward > 0.0)
                        {
                          cosmicval[count]=minval(diff_forward, diff_backward);
                          cosmicsig[count]=fabs(cosmicval[count])/limw[jjj*nx+iii];
                          ++count; //count number of cosmic rays (using first differences)
                          cmarr[jjj*nx+iii]=0;
                        }
                          if (diff_forward < 0.0 && diff_backward < 0.0)
                        {
                          cosmicval[count]=maxval(diff_forward, diff_backward);
                          cosmicsig[count]=fabs(cosmicval[count])/limw[jjj*nx+iii];
                          ++count; //count number of cosmic rays (using first differences)
                          cmarr[jjj*nx+iii]=0;
                        }

                          

                          

                          
                        }
                      }
                    }
            
                      }
                    }
                }
           
                }
              }
              }
    
          }
        
          
          //cosmic_ray_check[index_current]=count;
        
    




      //write out cosmic ray series

        //  if (cosmic_flag && kkk == 0){
        //    if (present[kkk] && keyvalue_fsn[kkk] >= fsn_first && keyvalue_fsn[kkk] <= fsn_last && time_fl[kkk] >= tfirst && time_fl[kkk] <= tlast)
        //      {
        //        printf("anormal writeout of first record\n");
        //        recout = dataout->records[kkk];
        //
        //      status=drms_setkey_int(recout, keyfsn, keyvalue_fsn[kkk]);
        //      status=drms_setkey_time(recout, keytobs, time_fl[kkk]);
        //      status=drms_setkey_int(recout, keycount, -1);
           //       status=drms_setkey_int(recout, fidkey, fid);
        //      status=drms_setkey_int(recout, keycamera, keyvalue_cam[kkk]);
        //      status=drms_setkey_int(recout, keyexmax, 0);
        //  status=drms_setkey_float(recout, keylimit, factor[cam][0]);
    
        //  drms_keyword_setdate(recout);
        
        //      if (keyvalue_cam[kkk] == cam_id_front) status=drms_setkey_string(recout, keyinstrument, camera_str_front);
        //  if (keyvalue_cam[kkk] == cam_id_side) status=drms_setkey_string(recout, keyinstrument, camera_str_side);
        //
        //  segout = drms_segment_lookup(recout, segmentname_cosmic);
            //  segout_val=drms_segment_lookup(recout, segmentname_val);
        //  segout_sig=drms_segment_lookup(recout, segmentname_sig);
        
     

        //  axisbad[0]=0;
        //  arrout=drms_array_create(type_int,1,axisbad,NULL,&status);
        //  cosmic_ray_data=arrout->data;
         
        //  arrout_val=drms_array_create(type_float,1,axisbad,NULL,&status);
        //  val_data=arrout_val->data;

        //  arrout_sig=drms_array_create(type_float,1,axisbad,NULL,&status);
        //  sig_data=arrout_sig->data;

        //  status=drms_segment_write(segout, arrout, 0);
        //  status=drms_segment_write(segout_val, arrout_val, 0);
        //  status=drms_segment_write(segout_sig, arrout_sig, 0);
       
        //  drms_free_array(arrout);
        //  drms_free_array(arrout_val);
        //  drms_free_array(arrout_sig);
        //   }
        //  }



        if (cosmic_flag && kkk >= 1){
          
         
          if (present[km1] && keyvalue_fsn[km1] >= fsn_first && keyvalue_fsn[km1] <= fsn_last && time_fl[km1] >= tfirst && time_fl[km1] <= tlast)
            {
              printf("fsn count %d %d\n", keyvalue_fsn[km1], count);
              printf("fsn order %d %d %d\n", fsnarr[last], fsnarr[current],fsnarr[next]);

              int cosfind, limit_flag;
              if (count == -1) cosfind=0; else if (count < limit_cosmic){cosfind=count; limit_flag=0;} else {cosfind=limit_cosmic; limit_flag=1;}


            recout = dataout->records[km1];

            status=drms_setkey_int(recout, keyfsn, keyvalue_fsn[km1]);
            status=drms_setkey_time(recout, keytobs, time_fl[km1]);
            status=drms_setkey_int(recout, keycount, count);
            status=drms_setkey_int(recout, fidkey, fid);
            status=drms_setkey_int(recout, keycamera, cam+1);
            status=drms_setkey_int(recout, keyexmax, limit_flag);
            status=drms_setkey_float(recout, keylimit, factor[cam][0]);
            drms_setkey_float (recout, "CRRMAX", crrmax);
    
            drms_keyword_setdate(recout);
        
            if (keyvalue_cam[km1] == cam_id_front) status=drms_setkey_string(recout, keyinstrument, camera_str_front);
            if (keyvalue_cam[km1] == cam_id_side) status=drms_setkey_string(recout, keyinstrument, camera_str_side);
        
            segout = drms_segment_lookup(recout, segmentname_cosmic);
            segout_val=drms_segment_lookup(recout, segmentname_val);
            segout_sig=drms_segment_lookup(recout, segmentname_sig);
        
     

            axisbad[0]=cosfind;
            arrout=drms_array_create(type_int,1,axisbad,NULL,&status);
            cosmic_ray_data=arrout->data;
         
            arrout_val=drms_array_create(type_float,1,axisbad,NULL,&status);
            val_data=arrout_val->data;

            arrout_sig=drms_array_create(type_float,1,axisbad,NULL,&status);
            sig_data=arrout_sig->data;


            for (i=0; i<cosfind; ++i){cosmic_ray_data[i]=cosmicarr[i]; val_data[i]=cosmicval[i]; sig_data[i]=cosmicsig[i];}

            status=drms_segment_write(segout, arrout, 0);
            status=drms_segment_write(segout_val, arrout_val, 0);
            status=drms_segment_write(segout_sig, arrout_sig, 0);
       
            drms_free_array(arrout);
            drms_free_array(arrout_val);
            drms_free_array(arrout_sig);
            }
        }

        
        

      



       
      /**************************************************************************************************/
      //permanent bad pixel calculation
        

        if (ccount >= cthreshold[cam] && kkk >= 1){
          update_stat=1;

          if (present_backward[km1] != 0) ++ccount1;
          if (present_forward[km1] != 0) ++ccount2;

#pragma omp parallel
          {
     
            if (present_backward[km1] != 0){
            
#pragma omp for private(jjj,iii)
              for (jjj=0; jjj<ny; ++jjj){
            for (iii=0; iii<nx; ++iii){
              //if (!isnan(arrimg[current][jjj][iii]) && arrimg[current][jjj][iii] > 0.0 && cmarr[jjj*nx+iii]){rhsp[jjj*nx+iii]=rhsp[jjj*nx+iii]+log((double)arrimg[current][jjj][iii]);  // //add up frames for backward
              if (!isnan(arrimg[current][jjj][iii]) && arrimg[current][jjj][iii] > 0.0){rhsp[jjj*nx+iii]=rhsp[jjj*nx+iii]+log((double)arrimg[current][jjj][iii]); 
                count_p[jjj*nx+iii]=count_p[jjj*nx+iii]+1;}
            }
              }
            }
    
     
            if (present_forward[km1] != 0){
            
        #pragma omp for private(jjj,iii)
              for (jjj=0; jjj<ny; ++jjj){
            for (iii=0; iii<nx; ++iii){
              //if (!isnan(arrimg[current][jjj][iii]) && arrimg[current][jjj][iii] > 0.0 && cmarr[jjj*nx+iii]){rhsm[jjj*nx+iii]=rhsm[jjj*nx+iii]+log((double)arrimg[current][jjj][iii]);   // //add up frames for forward
              if (!isnan(arrimg[current][jjj][iii]) && arrimg[current][jjj][iii] > 0.0){rhsm[jjj*nx+iii]=rhsm[jjj*nx+iii]+log((double)arrimg[current][jjj][iii]); 
                count_m[jjj*nx+iii]=count_m[jjj*nx+iii]+1;}
            }
              }
            }
          }

          
              
        }


        
        
          }

        /**************************************************/     
        //exchange indices

        holder=next;
        next=last;
        last=current;
        current=holder;
     


      }//end loop over filtergrams of same type
      
      /************************************************************************************/




      printtime();    
      printf("loop done\n");
      drms_close_records(dataout, DRMS_INSERT_RECORD);

      printf("reference flatfield and limb darkening\n");
   

      /*********************************/
      //calculate limb darkening function and correct for it

      limb_darkening((double)R_SUN, (double)XX0, (double)YY0, b_coef, b_order, limb_dark);

     
 #pragma omp parallel for private(jjj,iii)
      for (jjj=0; jjj<ny; ++jjj){
        for (iii=0; iii<nx; ++iii){
           
          if (count_p[jjj*nx+iii] > 0) rhsp[jjj*nx+iii]=(rhsp[jjj*nx+iii]/(double)count_p[jjj*nx+iii]-log(limb_dark[jjj*nx+iii])); 
          if (count_m[jjj*nx+iii] > 0) rhsm[jjj*nx+iii]=(rhsm[jjj*nx+iii]/(double)count_m[jjj*nx+iii]-log(limb_dark[jjj*nx+iii])); 
        }
      }

 
      //      if (debug)  //
      //{
      //printf("write debug\n");
      //printf("count_filtergram %d\n", count_filtergram);
          // FILE *fileptr;                                                                                                                          
          // int *aaa;                                                                                                                            
          // aaa=(int *)(malloc(count_filtergram*sizeof(int)));                                                                                              
          // fileptr = fopen ("/tmp20/richard/interpol/dd.bin", "w");                                                                                    
       // for (i=0;i<count_filtergram;i++){aaa[i]=farr[i];} fwrite ((char*)(aaa),sizeof(int),count_filtergram,fileptr);
       //  for (i=0;i<count_filtergram;i++){aaa[i]=barr[i];} fwrite ((char*)(aaa),sizeof(int),count_filtergram,fileptr);
          // fclose(fileptr);                                                                                                                            
      //    free(aaa);
      //}


      /*******************************************************************/
      //do flatfield calculations                                          /
      /*******************************************************************/

      printf("forward / backwards pairs %d %d \n", ccount1, ccount2);
 
      if (ccount1>= cthreshold[cam]){

  
      for (j=0; j<ny; ++j)  for (i=0; i<nx; ++i) flati[j*nx+i]=0.0;
 
      rotpairs=ccount1;

      cpa.norm=normconst[cam];
      cpa.omega=omegaconst[cam];
      cpa.convergence=convconst[cam];
      cpa.maxiter=maxiterconst[cam];

      if (flatfield_flag && status_flatfield_relative == 0)
    {
      status=flatfield(rhsp, rhsm, badpix, ccount, flati, param, cpa, deltat);
    }
      else
    {
      status=1;
    }
  

      if (status ==0)
    {
      /*********************************************************************/
      
      for (j=0; j<ny; ++j)  for (i=0; i<nx; ++i)flat[i][j]=exp(flati[j*nx+i]);
      update_stat=0;
   
      /*********************************************************************/
    }
      }
     


      if (status == 0) printf("flatfield for camera %d calculated \n", cam); else printf("could not calcuate flatfield for camera %d \n", cam);

    
     

      /*********************************************************************
      //write out  flatfield
      /********************************************************************/

 


      int axisout[2]={nx,ny};                             //size of the output array
      arrout_new  = drms_array_create(type_float,2,axisout,NULL,&status);
      printf("array size %ld\n", drms_array_size(arrout_new)/sizeof(float)); 
  
      flatfield_new=arrout_new->data;
  
 

      /*************************/
      //calculate apodization
             
       apod_circ(R_SUN*cpa.croprad*0.95, R_SUN*cpa.croprad*0.04, XX0-(float)nx/2.0, YY0-(float)ny/2.0, apod);

 
  
       /************************************/
       //copy and apodize flatfield

for (j=0; j<ny; ++j){
      for (i=0; i<nx; ++i){
    flatfield_new[j*nx+i]=(flat[i][j]-1.0)*apod[j*nx+i]+1.0;
      }
 }



 if (update_stat == 0)rot_new.rotpairs=rotpairs; else rot_new.rotpairs=0;
 rot_new.rotcadence=(float)deltat;
    //rot_new.flatfield_version=rot_cur.flatfield_version;


    
    printf("update_stat %d\n", update_stat);

    printf("query flat %s\n", query_flat);

    if (flatfield_flag)
      {
    status=write_flatfield_fid(arrout_new, cam+1, recnum_offpoint, t_0, focus, fid, rot_new);;
    if (status != 0){printf("could not write rotational flatfield\n"); return 1;}
      }


      drms_free_array(arrout_new);

    
    }
  else 
    {
      printf("No data records found\n");
      return 1;
    }     

    
  
  if (data != NULL)
    {
      printf("close DRMS session \n");
      drms_close_records(data,DRMS_FREE_RECORD);       
    }


  /************************************************************/
  //free arrays
  /****************************************************/


  free(param);
  free(rr);
  free(a1);
  free(a2);
  free(w1);
  free(w2);
  free(cmc);
  free(cosmicarr);
  free(cosmicsig);
  free(cosmicval);
  free(cmarr);
  free(apod);
      free(flat_off);
      free(flat_relative);
      free(badpix);
      free(app_flat);
      free(flati);

      free(rhsp);
      free(rhsm);
      free(count_p);
      free(count_m);
      free(limb_dark);

      free(limw);


  printtime();  
  printf("COMPLETED!\n");

  
  
  return 0;

}

/*
 *  Revision history
 *
 *  2017.10.31  Added optional argument for rmax_crfix to replace previously
 *  hard-coded value of 0.98
 *  Added setting of keyword CRRMAX if present for value of rmax_crfix
 * 
 */

//valgrind --leak-check=full --track-origins=yes --show-reachable=yes module_flatfield_256 input_series="su_production.lev1_hmi_test[986331-986430]" cosmic_rays=1

---