00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 #include <math.h>
00032
00033
00034
00035 #include <stdlib.h>
00036
00037 #include <cpl.h>
00038
00039 #include <irplib_utils.h>
00040
00041 #include <sinfo_tpl_utils.h>
00042 #include <sinfo_pfits.h>
00043 #include <sinfo_tpl_dfs.h>
00044 #include <sinfo_msg.h>
00045 #include <sinfo_error.h>
00046 #include <sinfo_utils_wrappers.h>
00047 #include <sinfo_globals.h>
00048 #include <sinfo_recipes.h>
00049 #include <sinfo_function_1d.h>
00050 #include <sinfo_functions.h>
00051 #include <sinfo_fit.h>
00052
00053
00054
00055
00056 static int sinfo_utl_table_test_create(cpl_plugin *) ;
00057 static int sinfo_utl_table_test_exec(cpl_plugin *) ;
00058 static int sinfo_utl_table_test_destroy(cpl_plugin *) ;
00059 static int sinfo_utl_table_test_shift(cpl_parameterlist *, cpl_frameset *) ;
00060
00061 static int
00062 sinfo_utl_table_test_amoeba_poly(cpl_parameterlist * parlist,
00063 cpl_frameset * framelist);
00064 static int
00065 sinfo_utl_table_test_amoeba_boltzmann(cpl_parameterlist * parlist,
00066 cpl_frameset * framelist);
00067
00068 static cpl_vector* sa_vx=NULL;
00069 static cpl_vector* sa_vy=NULL;
00070
00071 static double
00072 sinfo_fit_boltzmann(double p[]);
00073
00074
00075 static int
00076 sinfo_fitbkg(const double x[],
00077 const double a[],
00078 double *result);
00079
00080 static double
00081 sinfo_fac(const double x, const double t);
00082
00083 static double
00084 sinfo_fit_poly(double p[]);;
00085
00086 static cpl_table*
00087 sinfo_table_shift_column_spline3(cpl_table* t,
00088 const char* col,
00089 const double s);
00090
00091 static cpl_table*
00092 sinfo_table_shift_column_poly(cpl_table* t,
00093 const char* col,
00094 const double s,
00095 const int order);
00096
00097 static cpl_table*
00098 sinfo_table_shift_column_int(const cpl_table* t,
00099 const char* col,
00100 const double s,
00101 double* r);
00102
00103 static cpl_table*
00104 sinfo_table_shift_simple(cpl_table* inp,
00105 const char* col,
00106 const double shift);
00107
00108 #define NPOINT 1000
00109
00110
00111
00112
00113 static char sinfo_utl_table_test_description[] =
00114 "This recipe perform cubes combination.\n"
00115 "The input files are several cubeses\n"
00116 "their associated tags should be CUBE.\n"
00117 "The output is a cube PRO_CUBE resulting from the input cubes\n"
00118 "accurding to the value of op where op indicates the operation to be \n"
00119 "performed specified by the parameter sinfoni.sinfo_utl_table_test.op\n"
00120 "Information on relevant parameters can be found with\n"
00121 "esorex --params sinfo_utl_table_test\n"
00122 "esorex --help sinfo_utl_table_test\n"
00123 "\n";
00124
00125
00126
00127
00128
00132
00133
00135
00143
00144 int cpl_plugin_get_info(cpl_pluginlist * list)
00145 {
00146 cpl_recipe * recipe = cpl_calloc(1, sizeof *recipe ) ;
00147 cpl_plugin * plugin = &recipe->interface ;
00148
00149 cpl_plugin_init(plugin,
00150 CPL_PLUGIN_API,
00151 SINFONI_BINARY_VERSION,
00152 CPL_PLUGIN_TYPE_RECIPE,
00153 "sinfo_utl_table_test",
00154 "Combines a cube list in an output cube",
00155 sinfo_utl_table_test_description,
00156 "Andrea Modigliani",
00157 "Andrea.Modigliani@eso.org",
00158 sinfo_get_license(),
00159 sinfo_utl_table_test_create,
00160 sinfo_utl_table_test_exec,
00161 sinfo_utl_table_test_destroy) ;
00162
00163 cpl_pluginlist_append(list, plugin) ;
00164
00165 return 0;
00166 }
00167
00168
00177
00178 static int sinfo_utl_table_test_create(cpl_plugin * plugin)
00179 {
00180 cpl_recipe * recipe ;
00181 cpl_parameter * p ;
00182
00183
00184 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00185 recipe = (cpl_recipe *)plugin ;
00186 else return -1 ;
00187 cpl_error_reset();
00188 irplib_reset();
00189
00190
00191 recipe->parameters = cpl_parameterlist_new() ;
00192
00193
00194
00195 p = cpl_parameter_new_value("sinfoni.sinfo_utl_table_test.method",
00196 CPL_TYPE_STRING,
00197 "Output filename",
00198 "sinfoni.sinfo_utl_table_test",
00199 "sinfo_clean_mean");
00200 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "method") ;
00201 cpl_parameterlist_append(recipe->parameters, p) ;
00202
00203
00204 p = cpl_parameter_new_value("sinfoni.sinfo_utl_table_test.wshift",
00205 CPL_TYPE_DOUBLE,
00206 "Shift",
00207 "sinfoni.sinfo_utl_table_test",
00208 0.01) ;
00209
00210 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wshift") ;
00211 cpl_parameterlist_append(recipe->parameters, p) ;
00212
00213
00214 p = cpl_parameter_new_value("sinfoni.sinfo_utl_table_test.wstart",
00215 CPL_TYPE_DOUBLE,
00216 "Shift",
00217 "sinfoni.sinfo_utl_table_test",
00218 1.45) ;
00219
00220
00221 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wstart") ;
00222 cpl_parameterlist_append(recipe->parameters, p) ;
00223
00224
00225
00226 p = cpl_parameter_new_value("sinfoni.sinfo_utl_table_test.wend",
00227 CPL_TYPE_DOUBLE,
00228 "End Wavelength",
00229 "sinfoni.sinfo_utl_table_test",
00230 2.45) ;
00231
00232
00233 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wend") ;
00234 cpl_parameterlist_append(recipe->parameters, p) ;
00235
00236
00237
00238 return 0;
00239 }
00240
00241
00247
00248 static int sinfo_utl_table_test_exec(cpl_plugin * plugin)
00249 {
00250 cpl_recipe * recipe ;
00251 int status=0;
00252
00253 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00254 recipe = (cpl_recipe *)plugin ;
00255 else return -1 ;
00256 status+=sinfo_utl_table_test_shift(recipe->parameters, recipe->frames);
00257 status+=sinfo_utl_table_test_amoeba_poly(recipe->parameters, recipe->frames);
00258 status+=sinfo_utl_table_test_amoeba_boltzmann(recipe->parameters, recipe->frames);
00259
00260 return status ;
00261 }
00262
00263
00269
00270 static int sinfo_utl_table_test_destroy(cpl_plugin * plugin)
00271 {
00272 cpl_recipe * recipe ;
00273
00274
00275 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00276 recipe = (cpl_recipe *)plugin ;
00277 else return -1 ;
00278
00279 cpl_parameterlist_delete(recipe->parameters) ;
00280 return 0 ;
00281 }
00282
00283
00284
00291
00292 static int
00293 sinfo_utl_table_test_shift(cpl_parameterlist * parlist,
00294 cpl_frameset * framelist)
00295 {
00296
00297 cpl_table* t=NULL;
00298 cpl_table* t_shift=NULL;
00299 int np=NPOINT;
00300 cpl_parameter* p=NULL;
00301 const char* op=NULL;
00302 double ws=0;
00303 double we=0;
00304 double wd=0;
00305 double wshift=0;
00306 double pshift=0;
00307
00308 double* pw=NULL;
00309 double* pi=NULL;
00310
00311 int i=0;
00312
00313 double ra=0;
00314 const double pg=3.1415926535897932384626433832795;
00315
00316
00317 check(sinfo_dfs_set_groups(framelist),
00318 "Cannot identify RAW and CALIB frames") ;
00319
00320
00321 check_nomsg(p=cpl_parameterlist_find(parlist,
00322 "sinfoni.sinfo_utl_table_test.method"));
00323 check_nomsg(op=cpl_parameter_get_string(p));
00324
00325
00326 check_nomsg(p=cpl_parameterlist_find(parlist,
00327 "sinfoni.sinfo_utl_table_test.wstart"));
00328 check_nomsg(ws=cpl_parameter_get_double(p)) ;
00329
00330
00331 wd=(we-ws)/np;
00332 check_nomsg(p=cpl_parameterlist_find(parlist,
00333 "sinfoni.sinfo_utl_table_test.wend"));
00334 check_nomsg(we=cpl_parameter_get_double(p)) ;
00335
00336 check_nomsg(p=cpl_parameterlist_find(parlist,
00337 "sinfoni.sinfo_utl_table_test.wshift"));
00338 check_nomsg(wshift=cpl_parameter_get_double(p)) ;
00339
00340
00341
00342
00343 check_nomsg(t=cpl_table_new(np));
00344 check_nomsg(cpl_table_new_column(t,"WAVE",CPL_TYPE_DOUBLE));
00345 check_nomsg(cpl_table_new_column(t,"INT",CPL_TYPE_DOUBLE));
00346 check_nomsg(cpl_table_fill_column_window(t,"WAVE",0,np,0));
00347 check_nomsg(cpl_table_fill_column_window(t,"INT",0,np,0));
00348
00349 check_nomsg(pw=cpl_table_get_data_double(t,"WAVE"));
00350 check_nomsg(pi=cpl_table_get_data_double(t,"INT"));
00351
00352
00353
00354 for(i=0;i<np;i++) {
00355 pw[i]=ws+i*wd;
00356 ra=(double)i/(double)np*6.*pg;
00357 pi[i]=cos(ra);
00358 }
00359 pshift=wshift/wd;
00360 check_nomsg(cpl_table_save(t, NULL, NULL, "out_cosine.fits", 0));
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375 check_nomsg(t_shift=sinfo_table_shift_simple(t,"INT",pshift));
00376 check_nomsg(cpl_table_save(t_shift,NULL,NULL, "out_cosine_shift3.fits", 0));
00377 sinfo_free_table(&t);
00378
00379 return -1;
00380 cleanup:
00381
00382 sinfo_free_table(&t);
00383 return -1;
00384
00385
00386 }
00387
00388
00389
00390
00397
00398 static int
00399 sinfo_utl_table_test_amoeba_poly(cpl_parameterlist * parlist,
00400 cpl_frameset * framelist)
00401 {
00402
00403 cpl_table* t=NULL;
00404 int np=NPOINT;
00405 cpl_parameter* p=NULL;
00406 const char* op=NULL;
00407 double ws=0;
00408 double we=0;
00409 double wd=0;
00410 double wshift=0;
00411
00412 double* pw=NULL;
00413 double* pi=NULL;
00414 double* pf=NULL;
00415 int i=0;
00416 int j=0;
00417 double a[3];
00418 double p0[3];
00419
00420
00421
00422 const int MP=4;
00423 const int NP=3;
00424 double y[MP];
00425 double** ap=NULL;
00426 double FTOL=2e-6;
00427 int nfunc=0;
00428
00429 check(sinfo_dfs_set_groups(framelist),
00430 "Cannot identify RAW and CALIB frames") ;
00431
00432
00433 check_nomsg(p=cpl_parameterlist_find(parlist,
00434 "sinfoni.sinfo_utl_table_test.method"));
00435 check_nomsg(op=cpl_parameter_get_string(p));
00436
00437
00438 check_nomsg(p=cpl_parameterlist_find(parlist,
00439 "sinfoni.sinfo_utl_table_test.wstart"));
00440 check_nomsg(ws=cpl_parameter_get_double(p)) ;
00441
00442
00443 wd=(we-ws)/np;
00444 check_nomsg(p=cpl_parameterlist_find(parlist,
00445 "sinfoni.sinfo_utl_table_test.wend"));
00446 check_nomsg(we=cpl_parameter_get_double(p)) ;
00447
00448 check_nomsg(p=cpl_parameterlist_find(parlist,
00449 "sinfoni.sinfo_utl_table_test.wshift"));
00450 check_nomsg(wshift=cpl_parameter_get_double(p)) ;
00451
00452
00453
00454 check_nomsg(t=cpl_table_new(np));
00455 check_nomsg(cpl_table_new_column(t,"WAVE",CPL_TYPE_DOUBLE));
00456 check_nomsg(cpl_table_new_column(t,"INT",CPL_TYPE_DOUBLE));
00457 check_nomsg(cpl_table_fill_column_window(t,"WAVE",0,np,0));
00458 check_nomsg(cpl_table_fill_column_window(t,"INT",0,np,0));
00459
00460 check_nomsg(pw=cpl_table_get_data_double(t,"WAVE"));
00461 check_nomsg(pi=cpl_table_get_data_double(t,"INT"));
00462 ws=-2.;
00463 we=+4.;
00464 a[0]=+1.;
00465 a[1]=-2.;
00466 a[2]=1.;
00467 wd=(we-ws)/np;
00468 for(i=0;i<np;i++) {
00469 pw[i]=ws+i*wd;
00470
00471 pi[i]=a[0]+a[1]*pw[i]+a[2]*pw[i]*pw[i]+1.*rand()/RAND_MAX;
00472 }
00473
00474
00475 sa_vx=cpl_vector_wrap(np,cpl_table_get_data_double(t,"WAVE"));
00476 sa_vy=cpl_vector_wrap(np,cpl_table_get_data_double(t,"INT"));
00477
00478 ap=(double**) cpl_calloc(MP,sizeof(double*));
00479 for(i=0;i<MP;i++) {
00480 ap[i]=cpl_calloc(NP,sizeof(double));
00481 }
00482
00483 for(i=0;i<MP;i++) {
00484 for(j=0;j<NP;j++) {
00485 ap[i][j]=0;
00486 }
00487 }
00488
00489
00490
00491 ap[0][0]=-3; ap[0][1]=-3; ap[0][2]=-3;
00492 ap[1][0]=+3; ap[1][1]=-3; ap[1][2]=-3;
00493 ap[2][0]=+3; ap[2][1]=+3; ap[2][2]=-3;
00494 ap[3][0]=+3; ap[3][1]=-3; ap[3][2]=+3;
00495
00496 sinfo_msg("Before amoeba fit");
00497 for(i=0;i<MP;i++) {
00498 for(j=0;j<NP;j++) {
00499 sinfo_msg("ap[%d][%d]=%g",i,j,ap[i][j]);
00500 }
00501 }
00502
00503 for(i=0;i<MP;i++) {
00504 p0[0]=ap[i][0];
00505 p0[1]=ap[i][1];
00506 p0[2]=ap[i][2];
00507 y[i]=sinfo_fit_poly(p0);
00508 }
00509
00510 sinfo_fit_amoeba(ap,y,NP,FTOL,sinfo_fit_poly,&nfunc);
00511
00512 sinfo_msg("After amoeba fit");
00513 for(i=0;i<MP;i++) {
00514 for(j=0;j<NP;j++) {
00515 sinfo_msg("ap[%d][%d]=%g",i,j,ap[i][j]);
00516 }
00517 }
00518
00519 check_nomsg(cpl_table_new_column(t,"FIT",CPL_TYPE_DOUBLE));
00520 check_nomsg(cpl_table_fill_column_window(t,"FIT",0,np,0));
00521 check_nomsg(pf=cpl_table_get_data_double(t,"FIT"));
00522
00523 wd=(we-ws)/np;
00524 for(i=0;i<np;i++) {
00525 pw[i]=ws+i*wd;
00526 pf[i]=ap[0][0]+ap[0][1]*pw[i]+ap[0][2]*pw[i]*pw[i];
00527 }
00528
00529 check_nomsg(cpl_table_save(t,NULL,NULL, "out_amoeba_poly.fits", 0));
00530 cpl_vector_unwrap(sa_vx);
00531 cpl_vector_unwrap(sa_vy);
00532
00533 sinfo_free_table(&t);
00534 return -1;
00535 cleanup:
00536 cpl_vector_unwrap(sa_vx);
00537 cpl_vector_unwrap(sa_vy);
00538
00539 sinfo_free_table(&t);
00540 return -1;
00541
00542
00543 }
00544
00545
00552
00553 static int
00554 sinfo_utl_table_test_amoeba_boltzmann(cpl_parameterlist * parlist,
00555 cpl_frameset * framelist)
00556 {
00557
00558 cpl_table* t=NULL;
00559 int np=NPOINT;
00560 cpl_parameter* p=NULL;
00561 const char* op=NULL;
00562 double ws=0;
00563 double we=0;
00564 double wd=0;
00565 double wshift=0;
00566
00567 double* pw=NULL;
00568 double* pi=NULL;
00569 double* pf=NULL;
00570 int i=0;
00571 int j=0;
00572 double a[3];
00573 double p0[3];
00574
00575
00576
00577 const int MP=4;
00578 const int NP=3;
00579 double y[MP];
00580 double** ap=NULL;
00581 double FTOL=2e-6;
00582 int nfunc=0;
00583 double max=0;
00584
00585 double bkg_min=0;
00586 double bkg_max=0;
00587 double p0_min=0;
00588 double p0_max=0;
00589 double p1_min=0;
00590 double p1_max=0;
00591 double p2_min=0;
00592 double p2_max=0;
00593
00594
00595 check(sinfo_dfs_set_groups(framelist),
00596 "Cannot identify RAW and CALIB frames") ;
00597
00598
00599 check_nomsg(p=cpl_parameterlist_find(parlist,
00600 "sinfoni.sinfo_utl_table_test.method"));
00601 check_nomsg(op=cpl_parameter_get_string(p));
00602
00603
00604 check_nomsg(p=cpl_parameterlist_find(parlist,
00605 "sinfoni.sinfo_utl_table_test.wstart"));
00606 check_nomsg(ws=cpl_parameter_get_double(p)) ;
00607
00608
00609 wd=(we-ws)/np;
00610 check_nomsg(p=cpl_parameterlist_find(parlist,
00611 "sinfoni.sinfo_utl_table_test.wend"));
00612 check_nomsg(we=cpl_parameter_get_double(p)) ;
00613
00614 check_nomsg(p=cpl_parameterlist_find(parlist,
00615 "sinfoni.sinfo_utl_table_test.wshift"));
00616 check_nomsg(wshift=cpl_parameter_get_double(p)) ;
00617
00618
00619
00620 check_nomsg(t=cpl_table_new(np));
00621 check_nomsg(cpl_table_new_column(t,"WAVE",CPL_TYPE_DOUBLE));
00622 check_nomsg(cpl_table_new_column(t,"INT",CPL_TYPE_DOUBLE));
00623 check_nomsg(cpl_table_fill_column_window(t,"WAVE",0,np,0));
00624 check_nomsg(cpl_table_fill_column_window(t,"INT",0,np,0));
00625
00626 check_nomsg(pw=cpl_table_get_data_double(t,"WAVE"));
00627 check_nomsg(pi=cpl_table_get_data_double(t,"INT"));
00628 ws=1.4;
00629 we=2.5;
00630 wd=(we-ws)/np;
00631 a[2]=280;
00632 a[1]=55.817665;
00633 a[0]=548.77802;
00634 for(i=0;i<np;i++) {
00635 pw[i]=ws+i*wd;
00636 pi[i]=sinfo_fac(pw[i],a[2])*(1.+0.1*rand()/RAND_MAX);
00637 }
00638 check_nomsg(max=cpl_table_get_column_max(t,"INT"));
00639 check_nomsg(cpl_table_duplicate_column(t,"THERMAL",t,"INT"));
00640 check_nomsg(cpl_table_divide_scalar(t,"THERMAL",max));
00641 check_nomsg(cpl_table_multiply_scalar(t,"THERMAL",a[1]));
00642 check_nomsg(cpl_table_add_scalar(t,"THERMAL",a[0]));
00643
00644 check_nomsg(sa_vx=cpl_vector_wrap(np,cpl_table_get_data_double(t,"WAVE")));
00645 check_nomsg(sa_vy=cpl_vector_wrap(np,cpl_table_get_data_double(t,"THERMAL")));
00646
00647 ap=(double**) cpl_calloc(MP,sizeof(double*));
00648 for(i=0;i<MP;i++) {
00649 ap[i]=cpl_calloc(NP,sizeof(double));
00650 }
00651
00652 for(i=0;i<MP;i++) {
00653 for(j=0;j<NP;j++) {
00654 ap[i][j]=0;
00655 }
00656 }
00657
00658
00659 bkg_min=cpl_table_get_column_min(t,"THERMAL");
00660 bkg_max=cpl_table_get_column_max(t,"THERMAL");
00661
00662 p0_min=bkg_min/2;
00663 p0_max=bkg_min*2;
00664 p1_min=bkg_min/2.;
00665 p1_max=bkg_max*2;
00666 p1_min=200;
00667 p2_max=300;
00668
00669 ap[0][0]=p0_min; ap[0][1]=p1_min; ap[0][2]=p2_min;
00670 ap[1][0]=p0_max; ap[1][1]=p1_min; ap[1][2]=p2_min;
00671 ap[2][0]=p0_min; ap[2][1]=p1_max; ap[2][2]=p2_min;
00672 ap[3][0]=p0_min; ap[3][1]=p1_min; ap[3][2]=p2_max;
00673
00674 sinfo_msg("Before amoeba fit");
00675 for(i=0;i<MP;i++) {
00676 for(j=0;j<NP;j++) {
00677 sinfo_msg("ap[%d][%d]=%g",i,j,ap[i][j]);
00678 }
00679 }
00680
00681 for(i=0;i<MP;i++) {
00682 p0[0]=ap[i][0];
00683 p0[1]=ap[i][1];
00684 p0[2]=ap[i][2];
00685 y[i]=sinfo_fit_boltzmann(p0);
00686 }
00687
00688 sinfo_fit_amoeba(ap,y,NP,FTOL,sinfo_fit_boltzmann,&nfunc);
00689
00690 sinfo_msg("After amoeba fit");
00691 for(i=0;i<MP;i++) {
00692 for(j=0;j<NP;j++) {
00693 sinfo_msg("ap[%d][%d]=%g",i,j,ap[i][j]);
00694 }
00695 }
00696
00697 check_nomsg(cpl_table_new_column(t,"FIT",CPL_TYPE_DOUBLE));
00698 check_nomsg(cpl_table_fill_column_window(t,"FIT",0,np,0));
00699 check_nomsg(pf=cpl_table_get_data_double(t,"FIT"));
00700
00701 wd=(we-ws)/np;
00702 for(i=0;i<np;i++) {
00703 pw[i]=ws+i*wd;
00704 pf[i]=sinfo_fac(pw[i],ap[0][2]);
00705 }
00706 check_nomsg(max=cpl_table_get_column_max(t,"FIT"));
00707 check_nomsg(cpl_table_divide_scalar(t,"FIT",max));
00708 check_nomsg(cpl_table_multiply_scalar(t,"FIT",ap[0][1]));
00709 check_nomsg(cpl_table_add_scalar(t,"FIT",ap[0][0]));
00710
00711 check_nomsg(cpl_table_save(t,NULL,NULL, "out_amoeba_boltzmann.fits", 0));
00712
00713 cpl_vector_unwrap(sa_vx);
00714 cpl_vector_unwrap(sa_vy);
00715 sinfo_free_table(&t);
00716 return -1;
00717 cleanup:
00718
00719 cpl_vector_unwrap(sa_vx);
00720 cpl_vector_unwrap(sa_vy);
00721 sinfo_free_table(&t);
00722 irplib_error_dump(CPL_MSG_ERROR,CPL_MSG_ERROR);
00723 return -1;
00724
00725
00726 }
00727
00728
00729 static cpl_table*
00730 sinfo_table_shift_column_int(const cpl_table* t,
00731 const char* col,
00732 const double s,
00733 double* r)
00734 {
00735 cpl_table* out=NULL;
00736 int is=(int)s;
00737 int nrow=0;
00738 int i=0;
00739
00740 double* pi=NULL;
00741 double* po=NULL;
00742
00743 cknull(t,"null input table");
00744 out=cpl_table_duplicate(t);
00745 *r=s-is;
00746 nrow=cpl_table_get_nrow(t);
00747 pi=cpl_table_get_data_double(t,col);
00748 po=cpl_table_get_data_double(out,col);
00749 sinfo_msg("shifting of %d pixels",is);
00750 if(is > 0 ) {
00751 for(i=0;i<nrow;i++) {
00752 if( ((i-is) >=0) && ((i-is) < nrow)) {
00753 po[i-is]=pi[i];
00754 }
00755 }
00756 } else {
00757 for(i=nrow-1;i>-1;i--) {
00758 if( ((i-is) >=0) && ((i-is) < nrow)) {
00759 po[i-is]=pi[i];
00760 }
00761 }
00762 }
00763 return out;
00764 cleanup:
00765 sinfo_free_table(&out);
00766 return NULL;
00767
00768 }
00769
00770
00771
00772 static cpl_table*
00773 sinfo_table_shift_column_poly(cpl_table* t,
00774 const char* col,
00775 const double shift,
00776 const int order)
00777 {
00778 cpl_table* out=NULL;
00779 int nrow=0;
00780 int i=0;
00781 int flag=0;
00782 int n_points=0;
00783 int firstpos=0;
00784 int z=0;
00785 float eval=0;
00786 float sum=0;
00787 float new_sum=0;
00788 float* pi=NULL;
00789 float* po=NULL;
00790 float* spec=NULL ;
00791 float* corrected_spec=NULL ;
00792 float* xnum=NULL ;
00793 float* tableptr=NULL;
00794
00795 cknull(t,"null input table");
00796 if ( order <= 0 ) {
00797 sinfo_msg_error("wrong order of interpolation polynom given!") ;
00798 goto cleanup;
00799 }
00800
00801 out=cpl_table_duplicate(t);
00802
00803 nrow=cpl_table_get_nrow(t);
00804 cpl_table_cast_column(t,col,"FINT",CPL_TYPE_FLOAT);
00805 cpl_table_cast_column(out,col,"FINT",CPL_TYPE_FLOAT);
00806 pi=cpl_table_get_data_float(t,"FINT");
00807 po=cpl_table_get_data_float(out,"FINT");
00808
00809 n_points = order + 1 ;
00810 if ( n_points % 2 == 0 ) {
00811 firstpos = (int)(n_points/2) - 1 ;
00812 } else {
00813 firstpos = (int)(n_points/2) ;
00814 }
00815 spec=cpl_calloc(nrow,sizeof(float)) ;
00816 corrected_spec=cpl_calloc(nrow,sizeof(float)) ;
00817 xnum=cpl_calloc(order+1,sizeof(float)) ;
00818
00819 for ( i = 0 ; i < n_points ; i++ ) {
00820 xnum[i] = i ;
00821 }
00822
00823
00824 for(i=0;i<nrow;i++) {
00825 corrected_spec[i] = 0. ;
00826 }
00827
00828 sum = 0. ;
00829 for ( z = 0 ; z < nrow ; z++ ) {
00830 spec[z] = pi[z] ;
00831 if (isnan(spec[z]) ) {
00832 spec[z] = 0. ;
00833
00834 for ( i = z - firstpos ; i < z-firstpos+n_points ; i++ ) {
00835 if ( i < 0 ) continue ;
00836 if ( i >= nrow) continue ;
00837 corrected_spec[i] = ZERO ;
00838 }
00839 }
00840 if ( z != 0 && z != nrow - 1 ) {
00841 sum += spec[z] ;
00842 }
00843 }
00844
00845 new_sum = 0. ;
00846 for ( z = 0 ; z < nrow ; z++ ) {
00847
00848
00849
00850
00851
00852
00853 if (isnan(corrected_spec[z])) continue ;
00854 if ( z - firstpos < 0 ) {
00855 tableptr = &spec[0] ;
00856 eval = shift + z ;
00857 } else if ( z - firstpos + n_points >= nrow ) {
00858 tableptr = &spec[nrow - n_points] ;
00859 eval = shift + z + n_points - nrow ;
00860 } else {
00861 tableptr = &spec[z-firstpos] ;
00862 eval = shift + firstpos ;
00863 }
00864
00865 flag=0;
00866 corrected_spec[z]=sinfo_new_nev_ille(xnum,tableptr,order,eval,&flag);
00867 if ( z != 0 && z != nrow - 1 ) {
00868 new_sum += corrected_spec[z] ;
00869 }
00870 }
00871
00872
00873 for (z = 0 ; z < nrow ; z++ ) {
00874 if ( new_sum == 0. ) {
00875 new_sum = 1. ;
00876 }
00877 if ( z == 0 ) {
00878 po[z] = ZERO ;
00879 } else if ( z == nrow - 1 ) {
00880 po[z] = ZERO ;
00881 } else if ( isnan(corrected_spec[z]) ) {
00882 po[z] = ZERO ;
00883 } else {
00884 corrected_spec[z] *= sum / new_sum ;
00885 po[z] = corrected_spec[z] ;
00886 }
00887 }
00888 check_nomsg(cpl_table_erase_column(t,"FINT"));
00889 check_nomsg(cpl_table_erase_column(out,col));
00890 check_nomsg(cpl_table_cast_column(out,"FINT",col,CPL_TYPE_DOUBLE));
00891 check_nomsg(cpl_table_erase_column(out,"FINT"));
00892
00893 sinfo_free_float(&spec) ;
00894 sinfo_free_float(&corrected_spec) ;
00895 sinfo_free_float(&xnum) ;
00896
00897 return out;
00898 cleanup:
00899
00900
00901 sinfo_free_float(&spec) ;
00902 sinfo_free_float(&corrected_spec) ;
00903 sinfo_free_float(&xnum) ;
00904 sinfo_free_table(&out);
00905 return NULL;
00906
00907
00908 }
00909
00910
00911
00912
00913 static cpl_table*
00914 sinfo_table_shift_column_spline3(cpl_table* t,
00915 const char* col,
00916 const double shift)
00917 {
00918 cpl_table* out=NULL;
00919 int nrow=0;
00920 int i=0;
00921 int z=0;
00922
00923 float sum=0;
00924 float new_sum=0;
00925
00926 float* pi=NULL;
00927 float* po=NULL;
00928 float* eval=NULL;
00929 float* xnum=NULL;
00930 float* spec=NULL;
00931 float* corrected_spec=NULL;
00932
00933 cknull(t,"null input table");
00934 out=cpl_table_duplicate(t);
00935
00936 nrow=cpl_table_get_nrow(t);
00937 check_nomsg(cpl_table_cast_column(t,col,"FINT",CPL_TYPE_FLOAT));
00938 check_nomsg(cpl_table_cast_column(out,col,"FINT",CPL_TYPE_FLOAT));
00939 pi=cpl_table_get_data_float(t,"FINT");
00940 po=cpl_table_get_data_float(out,"FINT");
00941
00942
00943
00944 xnum=cpl_calloc(nrow,sizeof(float)) ;
00945
00946 for ( i = 0 ; i < nrow ; i++ ) {
00947 xnum[i] = i ;
00948 }
00949
00950 spec=cpl_calloc(nrow,sizeof(float)) ;
00951 corrected_spec=cpl_calloc(nrow,sizeof(float)) ;
00952 eval=cpl_calloc(nrow,sizeof(float)) ;
00953
00954 sum = 0. ;
00955 for ( z = 0 ; z < nrow ; z++ ) {
00956 spec[z] = pi[z] ;
00957 if (isnan(spec[z]) ) {
00958 for ( i = z-1 ; i <= z+1 ; i++ ) {
00959 if ( i < 0 ) continue ;
00960 if ( i >= nrow) continue ;
00961 corrected_spec[i] = ZERO ;
00962 }
00963 spec[z] = 0. ;
00964 }
00965 sum += spec[z] ;
00966 eval[z] = (float)shift+(float)z ;
00967 }
00968
00969 if ( -1 == sinfo_function1d_natural_spline(xnum,spec, nrow,
00970 eval,corrected_spec, nrow))
00971 {
00972 sinfo_msg_error("error in spline interpolation!") ;
00973 goto cleanup;
00974 }
00975
00976 new_sum = 0. ;
00977 for ( z = 0 ; z < nrow ; z++ ) {
00978 if ( isnan(corrected_spec[z]) ) {
00979 continue ;
00980 }
00981 new_sum += corrected_spec[z] ;
00982 }
00983
00984 for ( z = 0 ; z < nrow ; z++ ) {
00985 if ( new_sum == 0. ) new_sum =1. ;
00986 {
00987 if ( isnan(corrected_spec[z]) ) {
00988 po[z] = ZERO ;
00989 } else {
00990 corrected_spec[z] *= sum / new_sum ;
00991 po[z] = corrected_spec[z] ;
00992 }
00993 }
00994 }
00995
00996 sinfo_free_float(&xnum);
00997 sinfo_free_float(&spec) ;
00998 sinfo_free_float(&corrected_spec) ;
00999 sinfo_free_float(&eval) ;
01000
01001 check_nomsg(cpl_table_erase_column(t,"FINT"));
01002 check_nomsg(cpl_table_erase_column(out,col));
01003 check_nomsg(cpl_table_cast_column(out,"FINT",col,CPL_TYPE_DOUBLE));
01004 check_nomsg(cpl_table_erase_column(out,"FINT"));
01005
01006 return out;
01007 cleanup:
01008
01009 sinfo_free_float(&xnum);
01010 sinfo_free_float(&spec) ;
01011 sinfo_free_float(&corrected_spec) ;
01012 sinfo_free_float(&eval) ;
01013 sinfo_free_table(&out);
01014 return NULL;
01015
01016
01017 }
01018
01019
01020 static cpl_table*
01021 sinfo_table_shift_simple(cpl_table* inp,
01022 const char* col,
01023 const double shift)
01024 {
01025
01026 int nrow=0;
01027 cpl_table* out=NULL;
01028 int is=(int)shift;
01029 double ds=shift-is;
01030 double* pi=NULL;
01031 double* po=NULL;
01032 double m=0;
01033 int i=0;
01034 cknull(inp,"null input table");
01035
01036 check_nomsg(nrow=cpl_table_get_nrow(inp));
01037 check_nomsg(out=cpl_table_duplicate(inp));
01038 check_nomsg(cpl_table_fill_column_window(out,col,0,nrow,0));
01039 check_nomsg(pi=cpl_table_get_data_double(inp,col));
01040 check_nomsg(po=cpl_table_get_data_double(out,col));
01041
01042
01043 for(i=0;i<nrow;i++) {
01044 if((i+is)>0 && (i+is+1) < nrow) {
01045 m=pi[i+is+1]-pi[i+is];
01046 po[i]=pi[i+is]+m*ds;
01047 }
01048 }
01049 return out;
01050 cleanup:
01051 sinfo_free_table(&out);
01052 return NULL;
01053
01054 }
01055
01056
01057
01064
01065
01066 static double
01067 sinfo_fit_poly(double p[])
01068
01069 {
01070
01071 double* px=NULL;
01072 double* py=NULL;
01073
01074 int i=0;
01075 int np=0;
01076
01077 double fy=0;
01078 double chi2=0;
01079
01080 check_nomsg(px= cpl_vector_get_data(sa_vx));
01081 check_nomsg(py= cpl_vector_get_data(sa_vy));
01082 check_nomsg(np= cpl_vector_get_size(sa_vx));
01083
01084 for(i=0;i<np;i++) {
01085
01086 fy=p[0]+p[1]*px[i]+p[2]*px[i]*px[i];
01087 chi2+=(py[i]-fy)*(py[i]-fy);
01088 }
01089
01090 return chi2;
01091 cleanup:
01092 return -1;
01093
01094 }
01095
01096 static double
01097 sinfo_fit_boltzmann(double p[])
01098
01099 {
01100
01101 double* px=NULL;
01102 double* py=NULL;
01103 double* pv=NULL;
01104 cpl_vector* vtmp=NULL;
01105 double max=0;
01106 int i=0;
01107 int np=0;
01108
01109 double chi2=0;
01110
01111 check_nomsg(px= cpl_vector_get_data(sa_vx));
01112 check_nomsg(py= cpl_vector_get_data(sa_vy));
01113 check_nomsg(np= cpl_vector_get_size(sa_vx));
01114 check_nomsg(vtmp=cpl_vector_duplicate(sa_vy));
01115 check_nomsg(pv=cpl_vector_get_data(vtmp));
01116
01117 for(i=0;i<np;i++) {
01118 pv[i]=sinfo_fac(px[i],p[2]);
01119
01120 }
01121 check_nomsg(max=cpl_vector_get_max(vtmp));
01122 if(max> 0) {
01123 check_nomsg(cpl_vector_divide_scalar(vtmp,max));
01124 check_nomsg(cpl_vector_multiply_scalar(vtmp,p[1]));
01125 check_nomsg(cpl_vector_add_scalar(vtmp,p[0]));
01126 }
01127
01128
01129 for(i=0;i<np;i++) {
01130 chi2+=(py[i]-pv[i])*(py[i]-pv[i]);
01131 }
01132
01133 return chi2;
01134 cleanup:
01135 return -1;
01136
01137 }
01138
01139 static int
01140 sinfo_fitbkg(const double x[],
01141 const double a[],
01142 double *result)
01143 {
01144
01145 double fac = sinfo_fac(x[0],a[2]);
01146 *result = a[0]+a[1]*fac;
01147
01148 return 0;
01149 }
01150
01151 static double
01152 sinfo_fac(const double x, const double t)
01153 {
01154
01155 double c=14387.7;
01156
01157 return pow(x,-5.)/(exp(c/(x*fabs(t)))-1.);
01158 }
01159
01160