sinfo_utl_table_test.c

00001 /* $Id: sinfo_utl_table_test.c,v 1.2 2006/12/04 13:40:23 amodigli Exp $
00002  *
00003  * This file is part of the SINFONI Pipeline
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: amodigli $
00023  * $Date: 2006/12/04 13:40:23 $
00024  * $Revision: 1.2 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 #include <math.h>
00032 /*-----------------------------------------------------------------------------
00033                                 Includes
00034  ----------------------------------------------------------------------------*/
00035 #include <stdlib.h>
00036 /* cpl */
00037 #include <cpl.h>
00038 /* irplib */
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                             Functions prototypes
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                             Static variables
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                                 Functions code
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     /* Get the recipe out of the plugin */
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     /* Create the parameters list in the cpl_recipe object */
00191     recipe->parameters = cpl_parameterlist_new() ; 
00192 
00193     /* Fill the parameters list */
00194     /* --stropt */
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     /* --doubleopt */
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     /* --doubleopt */
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     /* --doubleopt */
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     /* Return */
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     /* Get the recipe out of the plugin */
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     /* Get the recipe out of the plugin */
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   // read input parameters
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   // create the table
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   sinfo_msg("shift1");
00363   cknull_nomsg(t_shift=sinfo_table_shift_column_int(t, "INT", pshift,&wrest));
00364   check_nomsg(cpl_table_save(t_shift,NULL,NULL,"out_cosine_shift1.fits", 0));
00365   sinfo_msg("shift2");
00366   sinfo_free_table(&t_shift);
00367   cknull_nomsg(t_shift=sinfo_table_shift_column_poly(t,"INT", pshift,2));
00368   check_nomsg(cpl_table_save(t_shift,NULL,NULL,"out_cosine_shift2.fits", 0));
00369   sinfo_msg("shift3");
00370   sinfo_free_table(&t_shift);
00371   cknull_nomsg(t_shift=sinfo_table_shift_column_spline3(t_shift,"INT",pshift));
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   // Amoeba fit:
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   // read input parameters
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   //Prepare a template function (polynomial)
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     //pi[i]=(pw[i]-a[0])*(pw[i]-a[0])+1.*rand()/RAND_MAX;
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   // Amoeba part
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   //ap[0][0]=-5;
00489   //ap[1][0]=-2;
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   // Amoeba fit:
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   // read input parameters
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   //Prepare a template function (polynomial)
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   // Amoeba part
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   //ap[0][0]=-5;
00658   //ap[1][0]=-2;
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   /* fill the xa[] array for the polint function */
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      * now determine the arrays of size n_points with which the
00849      * polynom is determined and determine the position eval
00850      * where the polynom is evaluated in polynomial interpolation.
00851      * Take care of the points near the row edges!
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   /* fill the output spectrum */
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   /* fill the xa[] array for the spline function */
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   /* now we do the spline interpolation*/
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   /* fill output imagelist */
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     //fy=(px[i]-p[0])*(px[i]-p[0]);
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     //sinfo_msg("x=%g p=%g",px[i],pv[i]);
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 

Generated on Wed Jan 17 08:33:45 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4