test_simulate.c

00001 /* $Id: test_simulate.c,v 1.31 2007/11/29 14:42:34 jmlarsen Exp $
00002  *
00003  * This file is part of the FORS Library
00004  * Copyright (C) 2002-2006 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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
00019  */
00020 
00021 /*
00022  * $Author: jmlarsen $
00023  * $Date: 2007/11/29 14:42:34 $
00024  * $Revision: 1.31 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <test_simulate.h>
00033 
00034 #include <fors_image.h>
00035 #include <fors_dfs.h>
00036 #include <fors_data.h>
00037 #include <fors_pfits.h>
00038 #include <fors_utils.h>
00039 
00040 #include <cpl.h>
00041 
00042 #include <math.h>
00043 
00054 static const int det_nx = 400;        /* Unbinned detector pixels */
00055 static const int det_ny = 400;
00056 static const int binx = 2;
00057 static const int biny = 2;
00058 static const double ron   = 4.0;       /* ADU */
00059 static const double conad = 0.78;   /* e- / ADU */
00060 
00061 static const double bias_avg = 200; /* ADU */
00062 static const double dark_avg = 50;  /* ADU */
00063 static const char *const instrume = "fors2";
00064 static const char *const chip_id = "Test chip 234";
00065 static const char *const read_clock = "200Kps/2ports/low_gain";
00066 
00074 static cpl_frame *
00075 frame_new(const char *filename, const char *tag, cpl_frame_group group)
00076 {
00077     cpl_frame *f = cpl_frame_new();
00078 
00079     cpl_frame_set_filename(f, filename);
00080     cpl_frame_set_tag     (f, tag);
00081     cpl_frame_set_group   (f, group);
00082 
00083     return f;
00084 }
00085 
00091 void
00092 create_standard_keys(cpl_propertylist *header, double exptime)
00093 {
00094     int nx = det_nx / binx;
00095     int ny = det_ny / biny;
00096 
00097     cpl_propertylist_update_string(header, "ESO DPR TYPE", "some");
00098     cpl_propertylist_update_string(header, "ESO TPL ID", "tpl id.");
00099     cpl_propertylist_update_string(header, "ESO INS COLL NAME", "collimator name");
00100     cpl_propertylist_update_string(header, "ARCFILE", "archive filename");
00101     
00102     cpl_propertylist_update_string(header, FORS_PFITS_INSTRUME, instrume);
00103     cpl_propertylist_update_string(header, FORS_PFITS_FILTER_NAME, "R_SPECIAL");
00104 
00105     cpl_propertylist_update_double(header, FORS_PFITS_AIRMASS_START, 1.156);    
00106     cpl_propertylist_update_double(header, FORS_PFITS_AIRMASS_END  , 1.619);    
00107     
00108     cpl_propertylist_update_int   (header, FORS_PFITS_DET_NX, nx);
00109     cpl_propertylist_update_int   (header, FORS_PFITS_DET_NY, ny);
00110     cpl_propertylist_update_int   (header, FORS_PFITS_BINX, binx);
00111     cpl_propertylist_update_int   (header, FORS_PFITS_BINY, biny);
00112     cpl_propertylist_update_int   (header, FORS_PFITS_OVERSCANX, 0);
00113     cpl_propertylist_update_int   (header, FORS_PFITS_OVERSCANY, 0);
00114     cpl_propertylist_update_int   (header, FORS_PFITS_PRESCANX, 0);
00115     cpl_propertylist_update_int   (header, FORS_PFITS_PRESCANY, 0);
00116 
00117     cpl_propertylist_update_double(header, FORS_PFITS_PIXSCALE, 0.126);
00118     
00119     cpl_propertylist_update_int(header, FORS_PFITS_OUTPUTS, 1);
00120     cpl_propertylist_update_double(header, FORS_PFITS_CONAD[0], conad);
00121     
00122     /* Convert RON to e- units */
00123     cpl_propertylist_update_double(header, FORS_PFITS_RON[0], ron*conad);
00124     cpl_propertylist_update_double(header, FORS_PFITS_EXPOSURE_TIME, exptime);
00125 
00126     cpl_propertylist_update_string(header, FORS_PFITS_CHIP_ID, chip_id);
00127     cpl_propertylist_update_string(header, FORS_PFITS_READ_CLOCK, read_clock);
00128     
00129     /* WCS info, fine tuned to match simulated catalogue */
00130     {
00131         struct {
00132             const char *name;
00133             double value;
00134         } data[] = {
00135             {"CRVAL1", 8.1368333333},
00136             {"CRVAL2", -46.9576388889},
00137             {"CRPIX1", 1},
00138             {"CRPIX2", 1},
00139 //            {"CD1_1", 5.81347849634012E-21}, 
00140 //            {"CD1_2", 9.49444444444444E-05},
00141 //            {"CD2_1",-9.49444444444444E-05}, 
00142 //            {"CD2_2",-5.81347849634012E-21},
00143             {"CD1_1", 5.81347849634012E-20}, 
00144             {"CD1_2", 9.49444444444444E-04},
00145             {"CD2_1",-9.49444444444444E-04}, 
00146             {"CD2_2",-5.81347849634012E-20},
00147             {"PV2_1", 1.0},
00148             {"PV2_2", 0.0},
00149             {"PV2_3", 42.0},
00150             {"PV2_4", 0.0}, 
00151             {"PV2_5", 0.0}
00152         };
00153         
00154         unsigned i;
00155         for (i = 0; i < sizeof(data) / sizeof(*data); i++) {
00156             cpl_propertylist_append_double(header, data[i].name, data[i].value);
00157         }
00158     }
00159     
00160     return;
00161 }
00162 
00163 
00164 #undef cleanup
00165 #define cleanup \
00166 do { \
00167     fors_image_delete(&bias); \
00168     cpl_propertylist_delete(header); \
00169 } while(0)
00170 
00180 cpl_frame *
00181 create_bias(const char *filename, const char *tag, cpl_frame_group group)
00182 {
00183     int nx = det_nx / binx;
00184     int ny = det_ny / biny;
00185     double exptime = 0.0;
00186     
00187     cpl_image *data     = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00188     cpl_image *variance = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00189     fors_image *bias = NULL;
00190     cpl_propertylist *header = cpl_propertylist_new();
00191     
00192     create_standard_keys(header, exptime);
00193     cpl_propertylist_erase(header, FORS_PFITS_FILTER_NAME); /* No filter */
00194 
00195     {
00196         int x, y;
00197         for (y = 1; y <= ny; y++)
00198             for (x = 1; x <= nx; x++) {
00199                 cpl_image_set(data   , x, y, 
00200                               (int)(bias_avg + ron*fors_rand_gauss() + 0.5));
00201                 cpl_image_set(variance, x, y, ron*ron);
00202             }
00203     }
00204     
00205     bias = fors_image_new(data, variance);
00206     fors_image_save(bias, header, filename);
00207 
00208     assure( !cpl_error_get_code(), return NULL,
00209             "Saving bias to %s failed", filename );
00210 
00211     cleanup;
00212     return frame_new(filename, tag, group);
00213 }
00214 
00215 #undef cleanup
00216 #define cleanup \
00217 do { \
00218     fors_image_delete(&dark); \
00219     cpl_propertylist_delete(header); \
00220 } while(0)
00221 
00232 cpl_frame *
00233 create_dark(const char *filename, const char *tag, cpl_frame_group group)
00234 {
00235     int nx = det_nx / binx;
00236     int ny = det_ny / biny;
00237     double exptime = 3.0;
00238     
00239     cpl_image *data     = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00240     cpl_image *variance = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00241     fors_image *dark = NULL;
00242     cpl_propertylist *header = cpl_propertylist_new();
00243     
00244     create_standard_keys(header, exptime);
00245 
00246     {
00247         int x, y;
00248         for (y = 1; y <= ny; y++)
00249             for (x = 1; x <= nx; x++) {
00250                 double var = ron*ron + dark_avg/conad;
00251                 
00252                 cpl_image_set(data   , x, y, 
00253                               (int)(bias_avg + dark_avg + 
00254                                     sqrt(var)*fors_rand_gauss() 
00255                                     + 0.5));
00256                 cpl_image_set(variance, x, y, var);
00257             }
00258     }
00259     
00260     dark = fors_image_new(data, variance);
00261     fors_image_save(dark, header, filename);
00262 
00263     assure( !cpl_error_get_code(), return NULL, 
00264             "Saving dark to %s failed", filename );
00265     
00266     cleanup;
00267     return frame_new(filename, tag, group);
00268 }
00269 
00270 
00271 #undef cleanup
00272 #define cleanup \
00273 do { \
00274     fors_image_delete(&sflat); \
00275     cpl_propertylist_delete(header); \
00276 } while(0)
00277 
00287 cpl_frame *
00288 create_screen_flat(const char *filename, const char *tag, cpl_frame_group group)
00289 {
00290     int nx = det_nx / binx;
00291     int ny = det_ny / biny;
00292     double exptime = 3.0;
00293 
00294     cpl_image *data     = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00295     cpl_image *variance = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00296     fors_image *sflat = NULL;
00297     cpl_propertylist *header = cpl_propertylist_new();
00298     
00299     create_standard_keys(header, exptime);
00300 
00301     {
00302         int x, y;
00303         for (y = 1; y <= ny; y++)
00304             for (x = 1; x <= nx; x++) {
00305                 double medium_scale_structure = 1000*(2 + sin(x*30.0/nx + y*30.0/ny));
00306                 double  large_scale_structure = 1000*(1 + sin(x*4.0/nx));
00307                 double flat = 5000 + medium_scale_structure + large_scale_structure;
00308                 double var = ron*ron + flat/conad;
00309                 
00310                 cpl_image_set(data   , x, y, 
00311                               (int)(bias_avg + flat + sqrt(var)*fors_rand_gauss() + 0.5));
00312                 cpl_image_set(variance, x, y, var);
00313             }
00314     }
00315     
00316     sflat = fors_image_new(data, variance);
00317     fors_image_save(sflat, header, filename);
00318 
00319     assure( !cpl_error_get_code(), return NULL, 
00320             "Saving screen flat to %s failed", filename );
00321     
00322     cleanup;
00323     return frame_new(filename, tag, group);
00324 }
00325 
00326 #undef cleanup
00327 #define cleanup \
00328 do { \
00329     fors_image_delete(&sflat); \
00330     cpl_propertylist_delete(header); \
00331 } while(0)
00332 
00343 cpl_frame *
00344 create_sky_flat(const char *filename, const char *tag, cpl_frame_group group,
00345                 double exptime)
00346 {
00347     int nx = det_nx / binx;
00348     int ny = det_ny / biny;
00349 
00350     cpl_image *data     = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00351     cpl_image *variance = cpl_image_new(nx, ny, FORS_IMAGE_TYPE);
00352     fors_image *sflat = NULL;
00353     cpl_propertylist *header = cpl_propertylist_new();
00354     
00355     create_standard_keys(header, exptime);
00356 
00357     {
00358         int x, y;
00359         for (y = 1; y <= ny; y++)
00360             for (x = 1; x <= nx; x++) {
00361                 double medium_scale_structure = 1000*(2 + sin(x*30.0/nx - y*10.0/ny));
00362                 double  large_scale_structure = 1000*(1 + sin(x*4.0/nx));
00363                 double flat = exptime*(5000 + 
00364                                        medium_scale_structure + 
00365                                        large_scale_structure);
00366                 double var = ron*ron + flat/conad;
00367                 
00368                 cpl_image_set(data   , x, y, 
00369                               (int)(bias_avg + flat + sqrt(var)*fors_rand_gauss() + 0.5));
00370                 cpl_image_set(variance, x, y, var);
00371             }
00372     }
00373     
00374     sflat = fors_image_new(data, variance);
00375     fors_image_save(sflat, header, filename);
00376 
00377     assure( !cpl_error_get_code(), return NULL, 
00378             "Saving sky flat to %s failed", filename );
00379 
00380     cleanup;
00381     return frame_new(filename, tag, group);
00382 }
00383 
00391 cpl_frame *
00392 create_standard(const char *filename, const char *tag, cpl_frame_group group)
00393 {
00394     // fixme: add stars
00395     double exptime = 1.0;
00396     return create_sky_flat(filename, tag, group, exptime);
00397 }
00398 
00399 
00400 #undef cleanup
00401 #define cleanup \
00402 do { \
00403     cpl_table_delete(t); \
00404 } while(0)
00405 
00413 cpl_frame *
00414 create_std_cat(const char *filename, const char *tag, cpl_frame_group group)
00415 {
00416     cpl_table *t;
00417     struct {
00418         double ra, dec;
00419         double magnitude, dmagnitude;
00420         double col;
00421         const char *name;
00422     } 
00423     data[] =  {
00424         {8.15958, -47.0347, 15.824000, 0.001, 0.8, "object 1"},
00425         {8.14792, -46.9664, 12.895000, 0.002, -0.2, ""},
00426         {8.15083, -47.0092, 12.861000, 0.003, -0.3, " dff bject 1"},
00427         {8.15583, -47.0222, 16.540001, 0.001, 0.7, "-9"},
00428         {8.17167, -47.10  , 11.970000, 0.005, 0.12, NULL},
00429         {8.14833, -47.0567, 13.861000, 0.003, -0.2, ""},
00430         {8.1475 , -47.0411, 13.903000, 0.001, -0.8, "dddddddobject 1"},
00431         {7.92542,  2.62917, 15.446000, 0.002, -0.2, "start 1"},
00432     };
00433 
00434     int N = sizeof(data) / sizeof(*data);
00435     int i;
00436 
00437     t = cpl_table_new(N);
00438     cpl_table_new_column(t, FORS_DATA_STD_RA , CPL_TYPE_DOUBLE);
00439     cpl_table_new_column(t, FORS_DATA_STD_DEC, CPL_TYPE_DOUBLE);
00440     cpl_table_new_column(t, FORS_DATA_STD_NAME, CPL_TYPE_STRING);
00441     
00442     for (i = 0; i < FORS_NUM_FILTER; i++) {
00443         if (!cpl_table_has_column(t, FORS_DATA_STD_MAG[i])) {
00444         cpl_table_new_column(t, FORS_DATA_STD_MAG[i], CPL_TYPE_FLOAT);
00445     }
00446         if (!cpl_table_has_column(t, FORS_DATA_STD_DMAG[i])) {
00447         cpl_table_new_column(t, FORS_DATA_STD_DMAG[i], CPL_TYPE_FLOAT);
00448     }
00449         if (!cpl_table_has_column(t, FORS_DATA_STD_COL[i])) {
00450             cpl_table_new_column(t, FORS_DATA_STD_COL[i], CPL_TYPE_FLOAT);
00451         }
00452     }
00453 
00454     for (i = 0; i < N; i++) {
00455         int j;
00456         
00457         cpl_table_set_double(t, FORS_DATA_STD_RA  , i, data[i].ra);
00458         cpl_table_set_double(t, FORS_DATA_STD_DEC , i, data[i].dec);
00459         cpl_table_set_string(t, FORS_DATA_STD_NAME, i, data[i].name);
00460         
00461         for (j = 0; j < FORS_NUM_FILTER; j++) {
00462             cpl_table_set_float (t, FORS_DATA_STD_MAG[j], i, data[i].magnitude);
00463             cpl_table_set_float (t, FORS_DATA_STD_DMAG[j], i, data[i].dmagnitude);
00464             cpl_table_set_float (t, FORS_DATA_STD_COL[j], i, data[i].col);
00465         }
00466     }
00467 
00468     cpl_table_save(t, NULL, NULL, filename, CPL_IO_DEFAULT);
00469     assure( !cpl_error_get_code(), return NULL, 
00470             "Failed to save standard catalogue to %s", filename );
00471     
00472     cleanup;
00473     return frame_new(filename, tag, group);
00474 }
00475 
00476 
00477 
00478 #undef cleanup
00479 #define cleanup \
00480 do { \
00481     cpl_table_delete(t); \
00482 } while(0)
00483 
00491 cpl_frame *
00492 create_phot_table(const char *filename, const char *tag, cpl_frame_group group)
00493 {
00494     cpl_table *t;
00495     struct {
00496         double ext_coeff, dext_coeff;
00497         double color_term, dcolor_term;
00498         double expected_zeropoint, dexpected_zeropoint;
00499     } 
00500     data[FORS_NUM_FILTER] = {
00501         {0.4  , 0.01, -0.076, 0.01, 20, 0.2},
00502         {0.05 , 0.01,  0.033, 0.01, 21.123456, 0.2},
00503         {0.1  , 0.01,  0.01 , 0.01, 22, 0.2},
00504         {0.09 , 0.01, -0.02 , 0.01, -18, 0.2},
00505         {0.2  , 0.01,  0.03 , 0.01, 0, 0.2},
00506         {0.000, 0.01, -0.04 , 0.01, 1.0, 0.2},
00507     };
00508     
00509     int N = sizeof(FORS_DATA_FILTERLIST) / sizeof(*FORS_DATA_FILTERLIST);
00510     int i;
00511     
00512     t = cpl_table_new(N);
00513     cpl_table_new_column(t, FORS_DATA_PHOT_FILTER   , CPL_TYPE_STRING);
00514     cpl_table_new_column(t, FORS_DATA_PHOT_EXTCOEFF , CPL_TYPE_DOUBLE);
00515     cpl_table_new_column(t, FORS_DATA_PHOT_DEXTCOEFF , CPL_TYPE_DOUBLE);
00516     cpl_table_new_column(t, FORS_DATA_PHOT_COLORTERM, CPL_TYPE_DOUBLE);
00517     cpl_table_new_column(t, FORS_DATA_PHOT_DCOLORTERM, CPL_TYPE_DOUBLE);
00518     cpl_table_new_column(t, FORS_DATA_PHOT_ZEROPOINT, CPL_TYPE_DOUBLE);
00519     cpl_table_new_column(t, FORS_DATA_PHOT_DZEROPOINT, CPL_TYPE_DOUBLE);
00520 
00521     /* For each filtername (e.g. "U_BESS") find the matching filter (e.g. FILTER_U) */
00522     for (i = 0; i < N; i++) {
00523         cpl_table_set_string(t, FORS_DATA_PHOT_FILTER   , i, FORS_DATA_FILTERLIST[i].name);
00524 
00525         unsigned j;
00526         for (j = 0; j < FORS_NUM_FILTER; j++) {
00527             if (FORS_DATA_FILTERLIST[i].filter == j) {
00528                 
00529                 cpl_table_set_double(t, FORS_DATA_PHOT_EXTCOEFF , i, data[j].ext_coeff);
00530                 cpl_table_set_double(t, FORS_DATA_PHOT_DEXTCOEFF , i, data[j].dext_coeff);
00531                 cpl_table_set_double(t, FORS_DATA_PHOT_COLORTERM, i, data[j].color_term);
00532                 cpl_table_set_double(t, FORS_DATA_PHOT_DCOLORTERM, i, data[j].dcolor_term);
00533                 cpl_table_set_double(t, FORS_DATA_PHOT_ZEROPOINT, i, data[j].expected_zeropoint);
00534                 cpl_table_set_double(t, FORS_DATA_PHOT_DZEROPOINT, i, data[j].dexpected_zeropoint);
00535             }
00536         }
00537     }
00538 
00539     cpl_table_save(t, NULL, NULL, filename, CPL_IO_DEFAULT);
00540     assure( !cpl_error_get_code(), return NULL, 
00541             "Failed to save photometry table to %s", filename );
00542     
00543     cleanup;
00544     return frame_new(filename, tag, group);
00545 }

Generated on Wed Sep 10 07:31:54 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6