JSOC_Documentation: proj/limbfit/apps/limbfit.h Source File

00001 /* I.Scholl "Wed Nov  6 09:42:51 HST 2013" 
00002 */
00003 
00004 #include <string.h>
00005 #include <stdio.h>
00006 #include <stdlib.h>
00007 #include <math.h>
00008 #include <unistd.h>
00009 #include <gsl/gsl_math.h>
00010 #include <gsl/gsl_deriv.h>
00011 #include <gsl/gsl_rng.h>                                      
00012 #include <gsl/gsl_randist.h>                                  
00013 #include <gsl/gsl_vector.h>                                   
00014 #include <gsl/gsl_blas.h>                                     
00015 #include <gsl/gsl_multifit_nlin.h>                            
00016 #include <gsl/gsl_errno.h> 
00017 #include <gsl/gsl_min.h>    
00018 #include <time.h>
00019 
00020 #include "fitsio.h"
00021 #include "nrutil.h"
00022 
00023 #include <jsoc_main.h>
00024 
00025 #include "expmax.h"
00026 #include "expfit.h"
00027 
00028 #define CODE_NAME       "limbfit"
00029 #define CODE_VERSION    "V4.0r9" 
00030 #define CODE_DATE       "Mon Aug 31 17:16:24 PDT 2015" 
00031 #define LOGMSG1         "LIMBFIT"
00032 #define JSD_NAME        "su_scholl.hmi_lf.jsd"
00033 
00034 //#define dsin  "hmi.lev1c_nrt[]"
00035 //#define dsout "su_scholl.limbfit"
00036 //#define LOG_DIR   "~/LOGS/"
00037 //#define TMP_DIR   "~/TMP/"
00038 
00039 #define NUMRECPERTRANS 72
00040 
00041 /*
00042 drms_open_records
00043 drms_create_records
00044 drms_array_create
00045 drms_segment_write
00046 drms_segment_write_from_file
00047 drms_set_key_string for the final status of the current processed record:
00048      (because if I can't write the final status of the processing even the record will be in a incoherent state...)
00049 */
00050 //---------------------------------------------------ERRORS
00051 //GENERAL FAILURES -> ABORT
00052 #define ERR_EXIT                            1
00053 #define ERR_USAGE                           -2
00054 #define ERR_MALLOC_FAILED                   -11
00055 #define ERR_SPECIAL                         -100
00056 #define ERR_DRMS_WRITE                      -200
00057 #define ERR_DRMS_READ                       -201
00058 #define WAR_DRMS_NORECORD                   201
00059 #define DEBUG_MSG                           999
00060 
00061 //GENERAL FAILURES -> write errors
00062 #define ERR_DRMS_WRITE_KW                   -300
00063 #define ERR_DRMS_READ_MISSING_DATA          -301
00064 #define ERR_DRMS_READ_MISSING_KW            -302
00065 #define ERR_DRMS_READ_MISSING_XYR_LF        -303
00066 #define ERR_FITSIO                          -350
00067 #define ERR_FITSIO2                         -351 // for mini FITS on error
00068 #define ERR_NR_STACK_TOO_SMALL              -352
00069 
00070 //LIMBFIT FAILED -> write errors
00071 #define ERR_LIMBFIT_FAILED                  -501
00072 #define ERR_LIMBFIT_FIT_FAILED              -502 // error on exit from Marcelo's fitting routines
00073 #define ERR_LIMBFIT_FLDF_FAILED             -503
00074 #define ERR_DISK_OUTSIDE_IMAGE              -511
00075 #define ERR_SIZE_ANN_TOO_BIG                -512
00076 #define ERR_GSL_FINMIN_SET_FAILED           -541
00077 #define ERR_GSL_FINMIN_NOMIN_FAILED         -542
00078 #define ERR_GSL_FINMIN_PRO_FAILED           -543
00079 #define ERR_GSL_GAUSSFIT_SET_FAILED         -551
00080 #define ERR_GSL_GAUSSFIT_FDFSOLVER_FAILED   -552
00081 
00082 //----------------------------- Processing status codes per record
00083 #define PROCSTAT_OK                     "OK"
00084 #define PROCSTAT_NOK                    "NOK"
00085 #define PROCSTAT_NO_LF_FAILED           "LF_FAILED"
00086 #define PROCSTAT_NO_LF_MISSVALS         "NO_LF_MISSVALS"
00087 #define PROCSTAT_NO_LF_DARKIMG          "NO_LF_DARKIMG"
00088 #define PROCSTAT_NO_LF_OPENLOOP         "NO_LF_OPENLOOP"
00089 #define PROCSTAT_NO_LF_DB_READ_PB       "NO_LF_DB_READ_PB"
00090 #define PROCSTAT_NO_LF_XYR_LF_MISSING   "NO_LF_XYR_LF_MISSING"
00091 #define PROCSTAT_NO_LF_DB_WRITE_PB      "NO_LF_DB_WRITE_PB"
00092 #define PROCSTAT_NO_LF_FITS_WRITE_PB    "NO_LF_FITS_WRITE_PB"
00093 
00094 //---------------------------------- LIMBFIT PARAMETERS
00095 #define ANNULUS_WIDTH 200                   // 
00096 #define MAX_SIZE_ANN_VARS 8000000           // ! must be the same value than JPT in fortran code !
00097 #define NUM_LDF 180                         // n/jang=NUM_LDF+1
00098 #define NUM_RADIAL_BINS 64                  // n/jprf
00099 #define NUM_AB_BINS 256                     // n/jreg
00100 #define LO_LIMIT 32.0                       // ! the sum of these 2 must be equal to ANNULUS_WIDTH 
00101 #define UP_LIMIT 32.0                       // 
00102 #define INC_X -4.0                          // 
00103 #define INC_Y -4.0                          // 
00104 #define NUM_FITPNTS 9                       // 2*NUM_FITPNTS<NUM_RADIAL_BINS
00105 #define GUESS_RANGE 8                       //
00106 #define NB_ITER 2                           //
00107 #define BAD_PIXEL_VALUE -2147483648.0
00108 //#define   SKIPGC 1                            // skip the guess estimation, use X0/YO_LF
00109 //#define   IFAC 0                              // skip the center calculation, use X0/YO_LF
00110 #define AHI 70000.0                         // 
00111 
00112 // alternate parameters for low LDF threshold - needed for roll analysis
00113 #define AHI2      30000.0
00114 #define LO_LIMIT2 24.0  
00115 #define UP_LIMIT2 24.0
00116 #define NB_ITER2  1
00117 
00118 //------------------------------------------------------
00119 
00120 typedef struct {
00121     float       *data;          // image to analyze
00122     int         img_sz0;
00123     int         img_sz1;
00124     int         cc;
00125     int         spe;
00126     int         iter;
00127     int         fldf;
00128     double      ix;
00129     double      iy;
00130     double      ir;
00131     //int       sav;
00132 } LIMBFIT_INPUT;
00133 
00134 typedef struct {
00135     float       *anls;          // annulus passed from one image to the next
00136     long        anls_nbpix;     // <=> jk
00137     float       *pf_anls;       //
00138     float       *pl_anls;       //
00139     int         is_firstobs;    // 0=yes, 1=no
00140 } LIMBFIT_IO_PUT;
00141 
00142 typedef struct {    // output files content
00143         
00144     // general keywords
00145     int         numext;
00146     float       cenx;
00147     float       ceny;
00148     double      radius;
00149     double      cmean;
00150     float       max_limb;
00151     int         quality;
00152     int         error1;
00153     int         error2;
00154         
00155     // result data
00156     float*      fits_ldfs_data;         // main table / segment
00157     float*      fits_fulldfs;           // extension #2
00158     float*      fits_alpha_beta;        // extension #0
00159     double*     fits_params;            // extension #1
00160 
00161     // info to describe extension dimensions
00162     long        fits_ldfs_naxis1;       //  ldf_nrow
00163     long        fits_ldfs_naxis2;       //  ldf_ncol
00164     long        fits_fldfs_nrows;       //  fldf_nrow
00165     long        fits_fldfs_tfields;     //  fldf_ncol
00166     long        fits_ab_nrows;          //  alpha_beta_nrow
00167     long        fits_ab_tfields;        //  alpha_beta_ncol
00168     long        fits_params_nrows;      //  params_nrow
00169     long        fits_params_tfields;    //  params_ncol
00170 
00171     // processing parameters to save
00172     int         ann_wd;
00173     long        mxszannv;
00174     int         nb_ldf;
00175     int         nb_rdb;
00176     int         nb_abb;
00177     double      up_limit;
00178     double      lo_limit;
00179     double      inc_x;
00180     double      inc_y;
00181     int         nfitpnts;
00182     int         nb_iter;
00183     int         cc;
00184     double      ahi;
00185     int         nb_fbins;
00186     int         fldfr;
00187     
00188     // extra for error management
00189     char*       dsin;
00190     char*       comment;
00191     char*       code_date;
00192     char*       code_version;
00193     char*       code_name;
00194     char*       bld_vers;
00195     
00196     // not to save
00197     FILE        *opf;
00198     char*       tmp_dir;
00199     char*       dsout;
00200     int         debug;
00201     
00202 } LIMBFIT_OUTPUT;
00203 
00204 // C functions
00205 void    close_on_error(DRMS_Record_t *record_in,DRMS_Record_t *record_out, DRMS_Array_t *data_array); //, FILE *opf);
00206 int     do_one_limbfit(unsigned int fsn, DRMS_Record_t *record_in,DRMS_Record_t *record_out, 
00207                     LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios, int *status);
00208 double  fin_min(double A[], double m, int range, int debug, FILE *fd);
00209 int     gaussfit(double y[], double t[],double sigma[], double A[], double erro[], long N, int degf, int debug, FILE *opf);
00210 void    get_sdate(char *sdate);
00211 int     get_set_kw(int typ, char *kw, char *kw_txt, unsigned int fsn, DRMS_Record_t *record_in,
00212                     DRMS_Record_t *record_out, fitsfile *outfptr, int tbf, LIMBFIT_OUTPUT *lfr, int *status);
00213 int     indexx(unsigned long n, float *arr, unsigned long *indx);
00214 void    lf_logmsg(char * type1, char * type2, int return_code, int status, char *message, char *code_name, FILE *opf);
00215 void    lf_logmsg4fitsio(char *log_msg,char *log_msg_code,char *kw,unsigned int fsn,int status, FILE *opf);
00216 int     limbfit(LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios);
00217 float   median(float * tmed, int siz);
00218 int     mk_fldfs(float cmx, float cmy, double radius, int naxis_row, int naxis_col, 
00219             long npixels, float *data, float **save_full_ldf, int *bins1, int *bins2, FILE *opf, int debug);
00220 int     process_n_records_fsn(char * open_dsname, LIMBFIT_INPUT *lfv, LIMBFIT_OUTPUT *lfr, LIMBFIT_IO_PUT *lfw, int *status);
00221 int     process_all_records_smpl(char * open_dsname, LIMBFIT_INPUT *lfv, LIMBFIT_OUTPUT *lfr, LIMBFIT_IO_PUT *lfw, int *status);
00222 int     write_mini_output(char * errcode, DRMS_Record_t *record_in,DRMS_Record_t *record_out,int tbf, LIMBFIT_OUTPUT *lfr);
00223 void    sav_b0(float *pf_sb0, float *pl_sb0, float *pf_b0);
00224 void    sum_b0(float *beta, float *pf_b0, float *pl_b0);
00225 int     sort(unsigned long n, float *arr);
00226 
00227 // fortran subroutine
00228 void    limb_(float *anls, long *jk, float *cmx, float *cmy, float *r, int *nitr, int *ncut,
00229                 float* rprf, float* lprf, float *rsi, float *rso, float *dx, float *dy, 
00230                 float* alph, float* beta, int *ifail, float* b0, int *centyp, float *lahi);