| 1 |
/* I.Scholl "Wed Nov 6 09:42:51 HST 2013" |
| 2 |
*/ |
| 3 |
|
| 4 |
#include <string.h> |
| 5 |
#include <stdio.h> |
| 6 |
#include <stdlib.h> |
| 7 |
#include <math.h> |
| 8 |
#include <unistd.h> |
| 9 |
#include <gsl/gsl_math.h> |
| 10 |
#include <gsl/gsl_deriv.h> |
| 11 |
#include <gsl/gsl_rng.h> |
| 12 |
#include <gsl/gsl_randist.h> |
| 13 |
#include <gsl/gsl_vector.h> |
| 14 |
#include <gsl/gsl_blas.h> |
| 15 |
#include <gsl/gsl_multifit_nlin.h> |
| 16 |
#include <gsl/gsl_errno.h> |
| 17 |
#include <gsl/gsl_min.h> |
| 18 |
#include <time.h> |
| 19 |
|
| 20 |
#include "fitsio.h" |
| 21 |
#include "nrutil.h" |
| 22 |
|
| 23 |
#include <jsoc_main.h> |
| 24 |
|
| 25 |
#include "expmax.h" |
| 26 |
#include "expfit.h" |
| 27 |
|
| 28 |
#define CODE_NAME "limbfit" |
| 29 |
#define CODE_VERSION "V4.0r9" |
| 30 |
#define CODE_DATE "Mon Aug 31 17:16:24 PDT 2015" |
| 31 |
#define LOGMSG1 "LIMBFIT" |
| 32 |
#define JSD_NAME "su_scholl.hmi_lf.jsd" |
| 33 |
|
| 34 |
//#define dsin "hmi.lev1c_nrt[]" |
| 35 |
//#define dsout "su_scholl.limbfit" |
| 36 |
//#define LOG_DIR "~/LOGS/" |
| 37 |
//#define TMP_DIR "~/TMP/" |
| 38 |
|
| 39 |
#define NUMRECPERTRANS 72 |
| 40 |
|
| 41 |
/* |
| 42 |
drms_open_records |
| 43 |
drms_create_records |
| 44 |
drms_array_create |
| 45 |
drms_segment_write |
| 46 |
drms_segment_write_from_file |
| 47 |
drms_set_key_string for the final status of the current processed record: |
| 48 |
(because if I can't write the final status of the processing even the record will be in a incoherent state...) |
| 49 |
*/ |
| 50 |
//---------------------------------------------------ERRORS |
| 51 |
//GENERAL FAILURES -> ABORT |
| 52 |
#define ERR_EXIT 1 |
| 53 |
#define ERR_USAGE -2 |
| 54 |
#define ERR_MALLOC_FAILED -11 |
| 55 |
#define ERR_SPECIAL -100 |
| 56 |
#define ERR_DRMS_WRITE -200 |
| 57 |
#define ERR_DRMS_READ -201 |
| 58 |
#define WAR_DRMS_NORECORD 201 |
| 59 |
#define DEBUG_MSG 999 |
| 60 |
|
| 61 |
//GENERAL FAILURES -> write errors |
| 62 |
#define ERR_DRMS_WRITE_KW -300 |
| 63 |
#define ERR_DRMS_READ_MISSING_DATA -301 |
| 64 |
#define ERR_DRMS_READ_MISSING_KW -302 |
| 65 |
#define ERR_DRMS_READ_MISSING_XYR_LF -303 |
| 66 |
#define ERR_FITSIO -350 |
| 67 |
#define ERR_FITSIO2 -351 // for mini FITS on error |
| 68 |
#define ERR_NR_STACK_TOO_SMALL -352 |
| 69 |
|
| 70 |
//LIMBFIT FAILED -> write errors |
| 71 |
#define ERR_LIMBFIT_FAILED -501 |
| 72 |
#define ERR_LIMBFIT_FIT_FAILED -502 // error on exit from Marcelo's fitting routines |
| 73 |
#define ERR_LIMBFIT_FLDF_FAILED -503 |
| 74 |
#define ERR_DISK_OUTSIDE_IMAGE -511 |
| 75 |
#define ERR_SIZE_ANN_TOO_BIG -512 |
| 76 |
#define ERR_GSL_FINMIN_SET_FAILED -541 |
| 77 |
#define ERR_GSL_FINMIN_NOMIN_FAILED -542 |
| 78 |
#define ERR_GSL_FINMIN_PRO_FAILED -543 |
| 79 |
#define ERR_GSL_GAUSSFIT_SET_FAILED -551 |
| 80 |
#define ERR_GSL_GAUSSFIT_FDFSOLVER_FAILED -552 |
| 81 |
|
| 82 |
//----------------------------- Processing status codes per record |
| 83 |
#define PROCSTAT_OK "OK" |
| 84 |
#define PROCSTAT_NOK "NOK" |
| 85 |
#define PROCSTAT_NO_LF_FAILED "LF_FAILED" |
| 86 |
#define PROCSTAT_NO_LF_MISSVALS "NO_LF_MISSVALS" |
| 87 |
#define PROCSTAT_NO_LF_DARKIMG "NO_LF_DARKIMG" |
| 88 |
#define PROCSTAT_NO_LF_OPENLOOP "NO_LF_OPENLOOP" |
| 89 |
#define PROCSTAT_NO_LF_DB_READ_PB "NO_LF_DB_READ_PB" |
| 90 |
#define PROCSTAT_NO_LF_XYR_LF_MISSING "NO_LF_XYR_LF_MISSING" |
| 91 |
#define PROCSTAT_NO_LF_DB_WRITE_PB "NO_LF_DB_WRITE_PB" |
| 92 |
#define PROCSTAT_NO_LF_FITS_WRITE_PB "NO_LF_FITS_WRITE_PB" |
| 93 |
|
| 94 |
//---------------------------------- LIMBFIT PARAMETERS |
| 95 |
#define ANNULUS_WIDTH 200 // |
| 96 |
#define MAX_SIZE_ANN_VARS 8000000 // ! must be the same value than JPT in fortran code ! |
| 97 |
#define NUM_LDF 180 // n/jang=NUM_LDF+1 |
| 98 |
#define NUM_RADIAL_BINS 64 // n/jprf |
| 99 |
#define NUM_AB_BINS 256 // n/jreg |
| 100 |
#define LO_LIMIT 32.0 // ! the sum of these 2 must be equal to ANNULUS_WIDTH |
| 101 |
#define UP_LIMIT 32.0 // |
| 102 |
#define INC_X -4.0 // |
| 103 |
#define INC_Y -4.0 // |
| 104 |
#define NUM_FITPNTS 9 // 2*NUM_FITPNTS<NUM_RADIAL_BINS |
| 105 |
#define GUESS_RANGE 8 // |
| 106 |
#define NB_ITER 2 // |
| 107 |
#define BAD_PIXEL_VALUE -2147483648.0 |
| 108 |
//#define SKIPGC 1 // skip the guess estimation, use X0/YO_LF |
| 109 |
//#define IFAC 0 // skip the center calculation, use X0/YO_LF |
| 110 |
#define AHI 70000.0 // |
| 111 |
|
| 112 |
// alternate parameters for low LDF threshold - needed for roll analysis |
| 113 |
#define AHI2 30000.0 |
| 114 |
#define LO_LIMIT2 24.0 |
| 115 |
#define UP_LIMIT2 24.0 |
| 116 |
#define NB_ITER2 1 |
| 117 |
|
| 118 |
//------------------------------------------------------ |
| 119 |
|
| 120 |
typedef struct { |
| 121 |
float *data; // image to analyze |
| 122 |
int img_sz0; |
| 123 |
int img_sz1; |
| 124 |
int cc; |
| 125 |
int spe; |
| 126 |
int iter; |
| 127 |
int fldf; |
| 128 |
double ix; |
| 129 |
double iy; |
| 130 |
double ir; |
| 131 |
//int sav; |
| 132 |
} LIMBFIT_INPUT; |
| 133 |
|
| 134 |
typedef struct { |
| 135 |
float *anls; // annulus passed from one image to the next |
| 136 |
long anls_nbpix; // <=> jk |
| 137 |
float *pf_anls; // |
| 138 |
float *pl_anls; // |
| 139 |
int is_firstobs; // 0=yes, 1=no |
| 140 |
} LIMBFIT_IO_PUT; |
| 141 |
|
| 142 |
typedef struct { // output files content |
| 143 |
|
| 144 |
// general keywords |
| 145 |
int numext; |
| 146 |
float cenx; |
| 147 |
float ceny; |
| 148 |
double radius; |
| 149 |
double cmean; |
| 150 |
float max_limb; |
| 151 |
int quality; |
| 152 |
int error1; |
| 153 |
int error2; |
| 154 |
|
| 155 |
// result data |
| 156 |
float* fits_ldfs_data; // main table / segment |
| 157 |
float* fits_fulldfs; // extension #2 |
| 158 |
float* fits_alpha_beta; // extension #0 |
| 159 |
double* fits_params; // extension #1 |
| 160 |
|
| 161 |
// info to describe extension dimensions |
| 162 |
long fits_ldfs_naxis1; // ldf_nrow |
| 163 |
long fits_ldfs_naxis2; // ldf_ncol |
| 164 |
long fits_fldfs_nrows; // fldf_nrow |
| 165 |
long fits_fldfs_tfields; // fldf_ncol |
| 166 |
long fits_ab_nrows; // alpha_beta_nrow |
| 167 |
long fits_ab_tfields; // alpha_beta_ncol |
| 168 |
long fits_params_nrows; // params_nrow |
| 169 |
long fits_params_tfields; // params_ncol |
| 170 |
|
| 171 |
// processing parameters to save |
| 172 |
int ann_wd; |
| 173 |
long mxszannv; |
| 174 |
int nb_ldf; |
| 175 |
int nb_rdb; |
| 176 |
int nb_abb; |
| 177 |
double up_limit; |
| 178 |
double lo_limit; |
| 179 |
double inc_x; |
| 180 |
double inc_y; |
| 181 |
int nfitpnts; |
| 182 |
int nb_iter; |
| 183 |
int cc; |
| 184 |
double ahi; |
| 185 |
int nb_fbins; |
| 186 |
int fldfr; |
| 187 |
|
| 188 |
// extra for error management |
| 189 |
char* dsin; |
| 190 |
char* comment; |
| 191 |
char* code_date; |
| 192 |
char* code_version; |
| 193 |
char* code_name; |
| 194 |
char* bld_vers; |
| 195 |
|
| 196 |
// not to save |
| 197 |
FILE *opf; |
| 198 |
char* tmp_dir; |
| 199 |
char* dsout; |
| 200 |
int debug; |
| 201 |
|
| 202 |
} LIMBFIT_OUTPUT; |
| 203 |
|
| 204 |
// C functions |
| 205 |
void close_on_error(DRMS_Record_t *record_in,DRMS_Record_t *record_out, DRMS_Array_t *data_array); //, FILE *opf); |
| 206 |
int do_one_limbfit(unsigned int fsn, DRMS_Record_t *record_in,DRMS_Record_t *record_out, |
| 207 |
LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios, int *status); |
| 208 |
double fin_min(double A[], double m, int range, int debug, FILE *fd); |
| 209 |
int gaussfit(double y[], double t[],double sigma[], double A[], double erro[], long N, int degf, int debug, FILE *opf); |
| 210 |
void get_sdate(char *sdate); |
| 211 |
int get_set_kw(int typ, char *kw, char *kw_txt, unsigned int fsn, DRMS_Record_t *record_in, |
| 212 |
DRMS_Record_t *record_out, fitsfile *outfptr, int tbf, LIMBFIT_OUTPUT *lfr, int *status); |
| 213 |
int indexx(unsigned long n, float *arr, unsigned long *indx); |
| 214 |
void lf_logmsg(char * type1, char * type2, int return_code, int status, char *message, char *code_name, FILE *opf); |
| 215 |
void lf_logmsg4fitsio(char *log_msg,char *log_msg_code,char *kw,unsigned int fsn,int status, FILE *opf); |
| 216 |
int limbfit(LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios); |
| 217 |
float median(float * tmed, int siz); |
| 218 |
int mk_fldfs(float cmx, float cmy, double radius, int naxis_row, int naxis_col, |
| 219 |
long npixels, float *data, float **save_full_ldf, int *bins1, int *bins2, FILE *opf, int debug); |
| 220 |
int process_n_records_fsn(char * open_dsname, LIMBFIT_INPUT *lfv, LIMBFIT_OUTPUT *lfr, LIMBFIT_IO_PUT *lfw, int *status); |
| 221 |
int process_all_records_smpl(char * open_dsname, LIMBFIT_INPUT *lfv, LIMBFIT_OUTPUT *lfr, LIMBFIT_IO_PUT *lfw, int *status); |
| 222 |
int write_mini_output(char * errcode, DRMS_Record_t *record_in,DRMS_Record_t *record_out,int tbf, LIMBFIT_OUTPUT *lfr); |
| 223 |
void sav_b0(float *pf_sb0, float *pl_sb0, float *pf_b0); |
| 224 |
void sum_b0(float *beta, float *pf_b0, float *pl_b0); |
| 225 |
int sort(unsigned long n, float *arr); |
| 226 |
|
| 227 |
// fortran subroutine |
| 228 |
void limb_(float *anls, long *jk, float *cmx, float *cmy, float *r, int *nitr, int *ncut, |
| 229 |
float* rprf, float* lprf, float *rsi, float *rso, float *dx, float *dy, |
| 230 |
float* alph, float* beta, int *ifail, float* b0, int *centyp, float *lahi); |
