proj/flatfield/apps/cosmic_ray_post.c

 /*
 * cosmic_ray_post - Post processing of cosmic ray series for false positives
 *
 */

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

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

#define kRecSetIn      "input_series" //name of the series containing the input filtergrams
#define datumn  "datum"
#define cameran "camera"
#define fsnname      "fsn"
#define fsnf_name "fsn_first"
#define fsnl_name "fsn_last"

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

ModuleArgs_t module_args[] =        
{
     {ARG_STRING, kRecSetIn, "",  "Input data series."},
  {ARG_STRING, "output_series", "hmi.cosmic_rays",  "Output data series"},
     {ARG_STRING, datumn, "yyyy.mm.dd", "Datum string"},
     {ARG_INT, "hour", "00", "Hour"},
     {ARG_INT, cameran, 0, "Camera"},
     {ARG_INT, fsnf_name, "0"},
     {ARG_INT, fsnl_name, "2147483647"},
     {ARG_END}
};



void printtime();
int is_element(int value, int *array, int n);

int DoIt(void)
{

#include "module_flatfield_const.h"

  /*********************/
  //read input parameters

  const char *inRecQuery = cmdparams_get_str(&cmdparams, kRecSetIn, NULL); //cmdparams is defined in jsoc_main.h
  const char *datum =  cmdparams_get_str(&cmdparams, datumn, NULL);
  const char *output_series = cmdparams_get_str (&cmdparams, "output_series", NULL);
 
 int cameraint = cmdparams_get_int(&cmdparams, cameran, NULL);
 int hour= cmdparams_get_int(&cmdparams, "hour", NULL);
 int fsn_first=cmdparams_get_int(&cmdparams, fsnf_name, NULL);
 int fsn_last=cmdparams_get_int(&cmdparams, fsnl_name, NULL);


 /**********************/
 //define variables


  int  status= DRMS_SUCCESS, stat=DRMS_SUCCESS;

  DRMS_Segment_t *seg, *seg_val, *seg_sig;
  DRMS_Segment_t *segout, *segout_val, *segout_sig;
  
  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;

  DRMS_RecordSet_t *data, *dataout;
  DRMS_Record_t  *recout, *rec0;
  DRMS_Array_t *arrout, *arrout_val, *arrout_sig;

  int *cosmic_ray_data;
  float *val_data, *sig_data;

  DRMS_Type_Value_t  keyvalue_time;

  int ct,ct_prev,ct_next,ctr;
  float *sig, *lev;
  int *hits;
  int axisbad[1];

  int i, j, k, c; //loop variables


  int axisout[2]={nx,ny};

  
  /**********************/
  //*CHECK WHETHER THE FLATFIELD OUTPUT SERIES EXISTS                                                                    */

    
       drms_series_exists(drms_env, output_series, &status);
      if (status == DRMS_ERROR_UNKNOWNSERIES)
    {
      printf("Output series %s doesn't exist\n", output_series);       //if the output series does not exit
      exit(EXIT_FAILURE);                                        //we exit the program
    } 
  

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

      /*****************************/
      //build query string

      char fnumb[2]={""};
      char ffnumb[2]={""};
      char timefirst[256]="";

      strcat(timefirst, datum);
      strcat(timefirst, "_");
      sprintf(fnumb, "%2.2d", hour);
      strcat(timefirst, fnumb);
      strcat(timefirst, ":00:00.00_TAI");
 
 

      char query0[256]="";
      TIME tfirst, tlast;

      printf("fsnfs %d %d\n", fsn_first, fsn_last);


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

 char timelast[256]="";
 strcat(timelast, datum);
 strcat(timelast, "_");
 sprintf(ffnumb, "%2.2d", hour+1);
 strcat(timelast, ffnumb);
 strcat(timelast, ":59:59.99_TAI");


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

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

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


 strcat(query0, inRecQuery);
 strcat(query0, "[");
 strcat(query0, datefirst);
 strcat(query0, "-");
 strcat(query0, datelast);
 strcat(query0, "][?CAMERA=");
 sprintf(fnumb, "%1.1d", cameraint);
 strcat(query0, fnumb);
 strcat(query0, "?]");

  printf("query string time: %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-2);
       strcat(query0, fsnf);
      strcat(query0, "-");
      char fsnl[10]={""};
      sprintf(fsnl, "%d", fsn_last+2);
       strcat(query0, fsnl);
       strcat(query0, "][?CAMERA=");
      char fnumb[2]={""};
     sprintf(fnumb, "%1.1d", cameraint);
     strcat(query0, fnumb);
      strcat(query0, "?]");
      printf("query string fsn: %s\n", query0);  
    }

  
  /****************************/
  //open records

 
  data     = drms_open_records(drms_env,query0,&stat);

  if (data == NULL){printf("can not open records\n"); exit(EXIT_FAILURE);}

  drms_stage_records(data, 1, 0);

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


  /******************************/
  //read out keywords


if (stat == DRMS_SUCCESS && data != NULL && nRecs > 0)
  {
  int keyvalue_fid[nRecs];
  int keyvalue_cam[nRecs];
  DRMS_Record_t *rec0[nRecs];
  int  *cosmic_rays[nRecs];
  int keyvalue_recnum[nRecs];
  DRMS_Array_t *arr_cosmic[nRecs], *arr_level[nRecs], *arr_sig[nRecs];
  short wave[nRecs];
  short pol[nRecs];
  int count[nRecs];
  double time_fl[nRecs];
  int keyvalue_fsn[nRecs];
  int keyvalue_flag[nRecs];
  float keyvalue_factor[nRecs];
  float *significance[nRecs], *level[nRecs];
  int new_count[nRecs];
  int npairs[nRecs];
  float crrmax[nRecs];


  printtime();
  printf("create_records ...");
  dataout = drms_create_records (drms_env, nRecs, output_series,
      DRMS_PERMANENT, &stat);
  printf("done\n");
  printtime();

  printf("begin keyword  reading loop\n");
for (k=0; k<nRecs; ++k)
  {
    rec0[k]=data->records[k];

   
   
    count[k]= drms_getkey_int(rec0[k],keycount,&status);
    keyvalue_fid[k] = drms_getkey_int(rec0[k],fidkey,&status);
    keyvalue_cam[k] = drms_getkey_int(rec0[k],keycamera,&status);
    keyvalue_fsn[k]= drms_getkey_int(rec0[k],keyfsn,&status);
    keyvalue_flag[k] = drms_getkey_int(rec0[k],keyexmax,&status);
    keyvalue_factor[k] = drms_getkey_float(rec0[k],keylimit,&status);

    keyvalue_time=drms_getkey(rec0[k], keytobs, &type_time, &status);   ;   //read "T_OBS" of filtergram
    time_fl[k] = keyvalue_time.time_val;

    wave[k]=((keyvalue_fid[k]-10000)/10-5)/2;
    pol[k]=(keyvalue_fid[k] % 10);    
    
    keyvalue_recnum[k]=drms_getkey_int(rec0[k], recnumkey,&status);
    crrmax[k] = drms_getkey_float (rec0[k], "CRRMAX", &status);

      }
 printtime();

 

 /***********************************/
 //data reading loop


printf("begin data  reading loop\n");

 for (k=0; k<nRecs; ++k)
   {

     seg = drms_segment_lookup(rec0[k], segmentname_cosmic);
     seg_val=drms_segment_lookup(rec0[k], segmentname_val);
     seg_sig=drms_segment_lookup(rec0[k], segmentname_sig);
        
     arr_cosmic[k]= drms_segment_read(seg, type_int, &status);
     arr_level[k]=drms_segment_read(seg_val, type_float, &status);
     arr_sig[k]=drms_segment_read(seg_sig, type_float, &status);


     cosmic_rays[k] = arr_cosmic[k]->data;
     significance[k]=arr_sig[k]->data;
     level[k]=arr_level[k]->data;

   }

 printtime();

 /******************************/
 //test loop

  int elem_prev, elem_next;

  int *cosmic_new;
  float *val_new;
  float *sig_new;

  int counter=0;
  int falsecounter=0, goodcounter=0;

  cosmic_new=(int *)(malloc(limit_cosmic*sizeof(int)));
  val_new=(float *)(malloc(limit_cosmic*sizeof(float)));
  sig_new=(float *)(malloc(limit_cosmic*sizeof(float)));




printf("begin test loop\n");

 for (k=0; k<nRecs; ++k)
   {
     

     if (time_fl[k] >= tfirst && time_fl[k] <= tlast  && keyvalue_fsn[k] >= fsn_first && keyvalue_fsn[k] <= fsn_last)
       {

     hits=cosmic_rays[k];
     sig=significance[k];
     lev=level[k];

     ct=minval(count[k], limit_cosmic);

     if (count[k] > 0 && k >=1 && k < (nRecs-1)) //at least one count, not first or last record -> do weeding
       {
         
         ct_prev=minval(count[k-1], limit_cosmic);
         ct_next=minval(count[k+1], limit_cosmic);
  

         ctr=0;

         for (c=0; c<ct; ++c)
           {
        
             ++counter;
         elem_prev=is_element(hits[c], cosmic_rays[k-1], ct_prev);
         elem_next=is_element(hits[c], cosmic_rays[k+1], ct_next);

         if (elem_prev == 0 && elem_next == 0)
           {
             
             cosmic_new[ctr]=hits[c];
             val_new[ctr]=lev[c];
             sig_new[ctr]=sig[c];
             ++ctr;
             ++goodcounter;
           }
         else
           {
             ++falsecounter;
           }
           }
                   
         new_count[k]=ctr;
       }
     else  //else, just copy
       {
     ctr=count[k];

     for (c=0; c<ct; ++c)
       {
                
         cosmic_new[c]=hits[c];
         val_new[c]=lev[c];
         sig_new[c]=sig[c];                                  
        
       }
           
     new_count[k]=ctr;
       }

     /*************************/
     //set keywords
     

     recout = dataout->records[k];

     status=drms_setkey_int(recout, keyfsn, keyvalue_fsn[k]);
     status=drms_setkey_time(recout, keytobs, time_fl[k]);
     
     status=drms_setkey_int(recout, fidkey, keyvalue_fid[k]);
     status=drms_setkey_int(recout, keycamera, keyvalue_cam[k]);
     status=drms_setkey_int(recout, keyexmax, keyvalue_flag[k]);
     status=drms_setkey_float(recout, keylimit, keyvalue_factor[k]); 
     drms_setkey_float (recout, "CRRMAX", crrmax[k]); 
    
     drms_keyword_setdate(recout);

     status=drms_setkey_int(recout, keycount, new_count[k]);
    
     if (keyvalue_cam[k] == 2) status=drms_setkey_string(recout, keyinstrument, camera_str_front);
     if (keyvalue_cam[k] == 1) status=drms_setkey_string(recout, keyinstrument, camera_str_side);
        
        
     
     /************************************/
     //write out data segments


     if (new_count[k] >= 0)
       {

     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]=new_count[k];
     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<new_count[k]; ++i){cosmic_ray_data[i]=cosmic_new[i]; val_data[i]=val_new[i]; sig_data[i]=sig_new[i];} //copy new data array

       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);
       }
       }
    
   }
 
 /****************/
 //free arrays

 for (k=0; k<nRecs; ++k)
   {
     drms_free_array(arr_level[k]);
     drms_free_array(arr_cosmic[k]);
     drms_free_array(arr_sig[k]);
   }

   
 drms_close_records(dataout, DRMS_INSERT_RECORD);


 drms_close_records(data,DRMS_FREE_RECORD);

 free(cosmic_new);
 free(sig_new);
 free(val_new);


 printf("correctly identified cosmic rays %d out of %d %f\n", goodcounter, (goodcounter+falsecounter), (float)goodcounter/(float)(falsecounter+goodcounter));
 
  }
 else 
    {
      printf("No data records found\n");
    }     

 
   



 printtime();

 printf("COMLETED!\n");
 return 0;

}





int is_element(int value, int *array, int n)
{

  int isel, ctr;

  ctr=0;
  isel=0;
  while (ctr < n && isel == 0)
    {
      isel=(value == array[ctr]);
      ++ctr;
    }

  return isel;
}


void printtime()    
// print time
{
  time_t timer, timerd;
  char *timestring;
  int i;

    timerd=time(&timer);
    timestring=ctime(&timer);
        for (i=0; i<24; ++i) printf("%c", *(timestring+i));
    printf("\n");
}

/*
 *  Revision history
 *
 *  2017.11.01  Added optional argument for output_series to replace previously
 *  hard-coded value of "hmi.cosmic_rays" (defined as filename_cosmic2_out)
 *  Added propagation of keyword CRRMAX if present from input_series
 * 
 */

---