Return to limbfit_tas.h CVS log

Up to [Development] / JSOC / proj / limbfit / apps

v6.01 jk correction

/* I.Scholl 
	#define CODE_VERSION 	"V6.1" 
	#define CODE_DATE 		"Fri Mar  6 22:30:42 PST 2015" 
*/

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_rng.h>                                      
#include <gsl/gsl_randist.h>                                  
#include <gsl/gsl_vector.h>                                   
#include <gsl/gsl_blas.h>                                     
#include <gsl/gsl_multifit_nlin.h>                            
#include <gsl/gsl_errno.h> 
#include <gsl/gsl_min.h>    
#include <time.h>

#include "nrutil.h"

#include <jsoc_main.h>

#include "expmax.h"
#include "expfit.h"

#define CODE_VERSION 	"V6r1" 
#define CODE_DATE 		"Fri Mar  6 22:30:42 PST 2015" 
#define LOGMSG1			"LIMBFITS"
#if defined(AB512) && AB512
#define CODE_NAME 		"limbfit_tas2(ab512)"
#define	JSD_NAME		"scholl_limbfit_tas_ab512.jsd"
#elif defined(AB1024) && AB1024
#define CODE_NAME 		"limbfit_tas2(ab1024)"
#define	JSD_NAME		"scholl_limbfit_tas_ab1024.jsd"
#else
#define CODE_NAME 		"limbfit_tas"
#define	JSD_NAME		"scholl_limbfit_tas.jsd"
#endif

//#define dsin	"hmi.lev1c_nrt[]"
//#define dsout	"su_scholl.limbfit"
//#define LOG_DIR	"~/LOGS/"
//#define TMP_DIR	"~/TMP/"

#define NUMRECPERTRANS 960 // must be equal to Unitsize in the JSD

/*
drms_open_records
drms_create_records
drms_array_create
drms_segment_write
drms_segment_write_from_file
drms_set_key_string for the final status of the current processed record (because if I can't write the final status of the processing even the record will be in a incoherent state...)
*/
//---------------------------------------------------ERRORS
//GENERAL FAILURES -> ABORT
#define ERR_EXIT 							1
#define ERR_USAGE 							-2
#define ERR_MALLOC_FAILED 					-11
#define ERR_SPECIAL							-100
#define ERR_DRMS_WRITE 						-200
#define ERR_DRMS_READ 						-201
#define ERR_DRMS_ARRAY_CREATE 				-202
#define ERR_DRMS_SEGMENT_WRITE 				-203
#define ERR_DRMS_SEGMENT_LOOKUPNUM			-204
#define WAR_DRMS_NORECORD 					201
#define WAR_DRMS_FREE_ARRAY					202
#define DEBUG_MSG 							999
#define	VOID								0

//GENERAL FAILURES -> write errors
#define ERR_DRMS_WRITE_KW 					-300
#define ERR_DRMS_READ_MISSING_DATA 			-301
#define ERR_DRMS_READ_MISSING_KW 			-302
#define ERR_DRMS_READ_MISSING_XYR_LF		-303
#define ERR_NR_STACK_TOO_SMALL 				-352

//LIMBFIT FAILED -> write errors
#define ERR_LIMBFIT_FAILED 					-501
#define ERR_LIMBFIT_FIT_FAILED 				-502 // error on exit from Marcelo's fitting routines
#define ERR_LIMBFIT_FLDF_FAILED				-503
#define ERR_DISK_OUTSIDE_IMAGE 				-511
#define ERR_SIZE_ANN_TOO_BIG 				-512
#define ERR_GSL_FINMIN_SET_FAILED 			-541
#define ERR_GSL_FINMIN_NOMIN_FAILED			-542
#define ERR_GSL_FINMIN_PRO_FAILED 			-543
#define ERR_GSL_GAUSSFIT_SET_FAILED 		-551
#define ERR_GSL_GAUSSFIT_FDFSOLVER_FAILED 	-552

//----------------------------- Processing status codes per record
#define PROCSTAT_OK						"OK"
#define PROCSTAT_NOK					"NOK"
#define PROCSTAT_NO_LF_FAILED			"LF_FAILED"
#define PROCSTAT_NO_LF_MISSVALS 		"NO_LF_MISSVALS"
#define PROCSTAT_NO_LF_DARKIMG 			"NO_LF_DARKIMG"
#define PROCSTAT_NO_LF_OPENLOOP 		"NO_LF_OPENLOOP"
#define PROCSTAT_NO_LF_DB_READ_PB 		"NO_LF_DB_READ_PB"
#define PROCSTAT_NO_LF_XYR_LF_MISSING 	"NO_LF_XYR_LF_MISSING"
#define PROCSTAT_NO_LF_DB_WRITE_PB 		"NO_LF_DB_WRITE_PB"

//---------------------------------- LIMBFIT PARAMETERS
#define ANNULUS_WIDTH 200					// 
#define MAX_SIZE_ANN_VARS 8000000			// ! must be the same value than JPT in fortran code !
#define NUM_LDF 180							// n/jang=NUM_LDF+1
#define NUM_RADIAL_BINS 64					// n/jprf
#if defined(AB512) && AB512
#define NUM_AB_BINS 512					// n/jreg
#elif defined(AB1024) && AB1024
#define NUM_AB_BINS 1024					// n/jreg
#else
#define NUM_AB_BINS 256					// n/jreg
#endif
#define LO_LIMIT 32.0						// ! the sum of these 2 must be equal to ANNULUS_WIDTH 
#define UP_LIMIT 32.0						// 
#define INC_X -4.0							// 
#define INC_Y -4.0							// 
#define NUM_FITPNTS 9						// 2*NUM_FITPNTS<NUM_RADIAL_BINS
#define GUESS_RANGE 8						//
#define NB_ITER 2							//
#define BAD_PIXEL_VALUE -2147483648.0
//#define	SKIPGC 1							// skip the guess estimation, use X0/YO_LF
//#define	IFAC 0								// skip the center calculation, use X0/YO_LF
#define	AHI 70000.0							// 

// alternate parameters for low LDF threshold - needed for roll analysis
#define	AHI2 	  30000.0
#define LO_LIMIT2 24.0	
#define UP_LIMIT2 24.0
#define NB_ITER2  1

//------------------------------------------------------

typedef struct {
	float		*data;			// image to analyze
	int			img_sz0;
	int			img_sz1;
	int			cc;
	int			spe;
	int			iter;
	int			fldf;
	double		ix;
	double		iy;
	double		ir;
	int			src;
	//int		sav;
	unsigned int fsn;
} LIMBFIT_INPUT;

typedef struct {
	float 		*anls;			// annulus passed from one image to the next
	long		anls_nbpix;		// <=> jk
	float		*pf_anls;		//
	float		*pl_anls;		//
	int			is_firstobs;	// 0=yes, 1=no
} LIMBFIT_IO_PUT;

typedef struct {	// output files content
		
	// general keywords
	int			numext;
	float		cenx;
	float		ceny;
	double		radius;
	double		cmean;
	float		max_limb;
	int			quality;
	int			error1;
	int			error2;
		
	// result data
	float*		fits_ldfs_data; 		// main table / segment
	float*		fits_fulldfs; 			// extension #2
	float*		fits_alpha_beta;  	 	// extension #0
	float		fits_as[6];
	float		fits_es[6];

	// info to describe extension dimensions
	int		fits_ldfs_naxis1;		//	ldf_nrow
	int		fits_ldfs_naxis2;		//	ldf_ncol
	int		fits_fldfs_nrows;		// 	fldf_nrow
	int 	fits_fldfs_tfields;		//	fldf_ncol
	int 	fits_ab_naxis1;			//	alpha_beta_nrow
	int		fits_ab_naxis2;			//	alpha_beta_ncol

	// processing parameters to save
	int			ann_wd;
	long		mxszannv;
	int			nb_ldf;
	int			nb_rdb;
	int			nb_abb;
	double		up_limit;
	double		lo_limit;
	double		inc_x;
	double		inc_y;
	int			nfitpnts;
	int			nb_iter;
	int			cc;
	double		ahi;
	int 		nb_fbins;
	int			fldfr;
	
	// extra for error management
	char*		dsin;
	char*		comment;
	char*		code_date;
	char*		code_version;
	char*		code_name;
	char*		bld_vers;
	char*		series_name;	
	
	// not to save
	FILE 		*opf;
	char*		tmp_dir;
	char*		dsout;
	int			debug;
	
} LIMBFIT_OUTPUT;

// C functions
void	close_on_error(DRMS_Record_t *record_in,DRMS_Record_t *record_out, DRMS_Array_t *data_array); //, FILE *opf);
int		do_one_limbfit(DRMS_Record_t *record_in,DRMS_Record_t *record_out, 
					LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios, int *status);
double	fin_min(double A[], double m, int range, int debug, FILE *fd);
int		gaussfit(double y[], double t[],double sigma[], double A[], double erro[], long N, int degf, int debug, FILE *opf);
void	get_sdate(char *sdate);
int		indexx(unsigned long n, float *arr, unsigned long *indx);
void	lf_logmsg(char * type1, char * type2, int return_code, int status, char *message, char *code_name, FILE *opf);
int		limbfit(LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios);
float	median(float * tmed, int siz);
int		mk_fldfs(float cmx, float cmy, double radius, int naxis_row, int naxis_col, 
			long npixels, float *data, float **save_full_ldf, int *bins1, int *bins2, FILE *opf, int debug);
int		process_n_records_fsn(char * open_dsname, LIMBFIT_INPUT *lfv, LIMBFIT_OUTPUT *lfr, LIMBFIT_IO_PUT *lfw, int *status);
int		process_all_records_smpl(char * open_dsname, LIMBFIT_INPUT *lfv, LIMBFIT_OUTPUT *lfr, LIMBFIT_IO_PUT *lfw, int *status);
int		write_mini_output(char * errcode, DRMS_Record_t *record_in,DRMS_Record_t *record_out,int tbf, LIMBFIT_OUTPUT *lfr);
int		write_lf_keywords(char * errcode, DRMS_Record_t *record_out, LIMBFIT_OUTPUT *results, int pass, int src);
void	sav_b0(float *pf_sb0, float *pl_sb0, float *pf_b0);
void	sum_b0(float *beta, float *pf_b0, float *pl_b0);
int		sort(unsigned long n, float *arr);

// fortran subroutine
void	limb_(float *anls, int *jk, float *cmx, float *cmy, float *r, int *nitr, int *ncut,
				float* rprf, float* lprf, float *rsi, float *rso, float *dx, float *dy, 
				float* alph, float* beta, int *ifail, float* b0, int *centyp, float *lahi);