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
00032
00033
00034
00035
00036 #include <string.h>
00037 #include <math.h>
00038
00039 #include <cpl.h>
00040
00041 #include "irplib_std.h"
00042 #include "irplib_irlist.h"
00043
00044
00048
00051
00052
00053
00054
00055 static int irplib_std_setactive(const char *) ;
00056 static const irplib_irstd * irplib_std_get_closest_star(double, double) ;
00057
00058
00059
00060
00061
00062
00073
00074 irplib_irstd * irplib_std_get_mag_one_cat(
00075 double ra,
00076 double dec,
00077 irplib_band band,
00078 const char * cat,
00079 double * mag)
00080 {
00081 const irplib_irstd * refstar ;
00082
00083
00084 if (cat == NULL) return NULL ;
00085 if (!strcmp(cat, "all")) return NULL ;
00086
00087
00088 irplib_std_setactive("none");
00089 irplib_std_setactive(cat) ;
00090 refstar = irplib_std_get_closest_star(ra, dec) ;
00091
00092
00093 if (refstar == NULL) return NULL ;
00094
00095
00096 switch (band) {
00097 case BAND_J:
00098 if (refstar->mag_J < 98.0) {
00099 *mag = (double)(refstar->mag_J) ;
00100 return (irplib_irstd*)refstar ;
00101 } else return NULL ;
00102 case BAND_H:
00103 if (refstar->mag_H < 98.0) {
00104 *mag = (double)(refstar->mag_H) ;
00105 return (irplib_irstd*)refstar ;
00106 } else return NULL ;
00107 case BAND_K:
00108 if (refstar->mag_K < 98.0) {
00109 *mag = (double)(refstar->mag_K) ;
00110 return (irplib_irstd*)refstar ;
00111 } else return NULL ;
00112 case BAND_KS:
00113 if (refstar->mag_Ks < 98.0) {
00114 *mag = (double)(refstar->mag_Ks) ;
00115 return (irplib_irstd*)refstar ;
00116 } else return NULL ;
00117 case BAND_L:
00118 if (refstar->mag_L < 98.0) {
00119 *mag = (double)(refstar->mag_L) ;
00120 return (irplib_irstd*)refstar ;
00121 } else return NULL ;
00122 case BAND_M:
00123 if (refstar->mag_M < 98.0) {
00124 *mag = (double)(refstar->mag_M) ;
00125 return (irplib_irstd*)refstar ;
00126 } else return NULL ;
00127 case BAND_LP:
00128 if (refstar->mag_Lp < 98.0) {
00129 *mag = (double)(refstar->mag_Lp) ;
00130 return (irplib_irstd*)refstar ;
00131 } else return NULL ;
00132 case BAND_MP:
00133 if (refstar->mag_Mp < 98.0) {
00134 *mag = (double)(refstar->mag_Mp) ;
00135 return (irplib_irstd*)refstar ;
00136 } else return NULL ;
00137 default:
00138 return NULL ;
00139 }
00140 }
00141
00142
00155
00156 cpl_vector * irplib_std_get_conversion(
00157 const cpl_bivector * spec,
00158 double dit,
00159 double surface,
00160 double gain,
00161 double mag)
00162 {
00163 double h = 6.62e-27 ;
00164 double c = 3e18 ;
00165 cpl_vector * wave ;
00166 cpl_vector * extr ;
00167 cpl_vector * out ;
00168 double factor ;
00169
00170
00171 if (spec == NULL) return NULL ;
00172 if (dit <= 0.0) return NULL ;
00173
00174
00175 wave = cpl_bivector_get_x(spec) ;
00176 extr = cpl_bivector_get_y(spec) ;
00177
00178
00179 out = cpl_vector_duplicate(extr) ;
00180
00181
00182 cpl_vector_divide_scalar(out, dit) ;
00183
00184
00185 cpl_vector_divide_scalar(out, surface) ;
00186
00187
00188 cpl_vector_multiply_scalar(out, gain) ;
00189
00190
00191 factor = pow(10, mag/2.5) ;
00192 cpl_vector_multiply_scalar(out, factor) ;
00193
00194
00195 factor = (cpl_vector_get(wave, cpl_vector_get_size(wave)-1) -
00196 cpl_vector_get(wave, 0)) / cpl_vector_get_size(wave) ;
00197 cpl_vector_divide_scalar(out, factor) ;
00198
00199
00200 cpl_vector_multiply_scalar(out, h*c) ;
00201 cpl_vector_divide(out, wave) ;
00202
00203 return out ;
00204 }
00205
00206
00214
00215 cpl_vector * irplib_std_get_mag_zero(
00216 const cpl_bivector * sed,
00217 const cpl_vector * waves,
00218 double cent_wl)
00219 {
00220 double wmin, wmax, wstep ;
00221 int nb_sed ;
00222 double * sed_x ;
00223 double * sed_y ;
00224 cpl_bivector * sed_loc ;
00225 double * sed_loc_x ;
00226 double * sed_loc_y ;
00227 cpl_vector * out ;
00228 cpl_bivector * out_biv ;
00229 double f0_jan, f0_erg, cent_val ;
00230 int i ;
00231
00232
00233 if (sed == NULL) return NULL ;
00234 if (waves == NULL) return NULL ;
00235
00236
00237 nb_sed = cpl_bivector_get_size(sed) ;
00238 sed_x = cpl_bivector_get_x_data(sed) ;
00239 sed_y = cpl_bivector_get_y_data(sed) ;
00240 wstep = sed_x[1] - sed_x[0] ;
00241 wmin = cpl_vector_get(waves, 0) ;
00242 wmax = cpl_vector_get(waves, cpl_vector_get_size(waves)-1) ;
00243
00244
00245 sed_loc = cpl_bivector_new(nb_sed + 4) ;
00246 sed_loc_x = cpl_bivector_get_x_data(sed_loc) ;
00247 sed_loc_y = cpl_bivector_get_y_data(sed_loc) ;
00248 for (i=0 ; i<nb_sed ; i++) {
00249 sed_loc_x[i+2] = sed_x[i] ;
00250 sed_loc_y[i+2] = sed_y[i] ;
00251 }
00252
00253
00254 sed_loc_x[1] = sed_loc_x[2] - wstep ;
00255 if (sed_loc_x[2] < wmin) {
00256 sed_loc_x[0] = sed_loc_x[1] - wstep ;
00257 } else {
00258 sed_loc_x[0] = wmin - wstep ;
00259 }
00260 sed_loc_y[0] = 1e20 ;
00261 sed_loc_y[1] = 1e20 ;
00262
00263
00264 sed_loc_x[nb_sed+2] = sed_loc_x[nb_sed+1] + wstep ;
00265 if (sed_loc_x[nb_sed+1] > wmax) {
00266 sed_loc_x[nb_sed+3] = sed_loc_x[nb_sed+2] + wstep ;
00267 } else {
00268 sed_loc_x[nb_sed+3] = wmax + wstep ;
00269 }
00270 sed_loc_y[nb_sed+2] = 1e20 ;
00271 sed_loc_y[nb_sed+3] = 1e20 ;
00272
00273
00274 out = cpl_vector_duplicate(waves) ;
00275 out_biv = cpl_bivector_wrap_vectors((cpl_vector*)waves, out) ;
00276
00277
00278 if (cpl_bivector_interpolate_linear(out_biv, sed_loc) != CPL_ERROR_NONE) {
00279 cpl_msg_error(cpl_func, "Cannot interpolate the wavelength") ;
00280 cpl_bivector_unwrap_vectors(out_biv) ;
00281 cpl_vector_delete(out) ;
00282 cpl_bivector_delete(sed_loc) ;
00283 return NULL ;
00284 }
00285 cpl_bivector_unwrap_vectors(out_biv) ;
00286 cpl_bivector_delete(sed_loc) ;
00287
00288
00289 f0_jan = 5513.15 / ( pow(cent_wl,3) * (exp(1.2848/cent_wl)-1) ) ;
00290
00291
00292 f0_erg = f0_jan * 1e-26 * 1e7 * 3e18 / (1e4 * cent_wl*cent_wl*1e4*1e4) ;
00293
00294
00295 cent_val = cpl_vector_get(out, cpl_vector_get_size(out)/2) ;
00296 if (cent_val <= 0.0) {
00297 cpl_msg_error(cpl_func, "Negative or 0 central value") ;
00298 cpl_vector_delete(out) ;
00299 return NULL ;
00300 }
00301 cpl_vector_multiply_scalar(out, f0_erg/cent_val) ;
00302
00303
00304 return out ;
00305 }
00306
00307
00317
00318 cpl_bivector * irplib_std_get_sed(
00319 const char * seds_file,
00320 const char * sptype)
00321 {
00322 cpl_table * seds ;
00323 cpl_bivector * out ;
00324 cpl_vector * wave ;
00325 cpl_vector * sed ;
00326 cpl_bivector * tmp ;
00327 int nlines ;
00328
00329
00330 if (seds_file == NULL) return NULL ;
00331 if (sptype == NULL) return NULL ;
00332
00333
00334 if ((seds = cpl_table_load(seds_file, 1, 0)) == NULL) {
00335 cpl_msg_error(cpl_func, "Cannot load the table") ;
00336 return NULL ;
00337 }
00338
00339
00340 if (!cpl_table_has_column(seds, sptype)) {
00341 cpl_msg_error(cpl_func, "SED of the requested star not available") ;
00342 cpl_table_delete(seds) ;
00343 return NULL ;
00344 }
00345
00346
00347 nlines = cpl_table_get_nrow(seds) ;
00348
00349
00350 if ((wave = cpl_vector_wrap(nlines,
00351 cpl_table_get_data_double(seds, "Wavelength"))) == NULL) {
00352 cpl_msg_error(cpl_func, "Cannot get the Wavelength column") ;
00353 cpl_table_delete(seds) ;
00354 return NULL ;
00355 }
00356
00357
00358 if ((sed = cpl_vector_wrap(nlines,
00359 cpl_table_get_data_double(seds, sptype))) == NULL) {
00360 cpl_msg_error(cpl_func, "Cannot get the SED column") ;
00361 cpl_table_delete(seds) ;
00362 cpl_vector_unwrap(wave) ;
00363 return NULL ;
00364 }
00365 tmp = cpl_bivector_wrap_vectors(wave, sed) ;
00366
00367
00368 out = cpl_bivector_duplicate(tmp) ;
00369
00370
00371 cpl_bivector_unwrap_vectors(tmp) ;
00372 cpl_vector_unwrap(wave) ;
00373 cpl_vector_unwrap(sed) ;
00374 cpl_table_delete(seds) ;
00375
00376
00377 return out ;
00378 }
00379
00380
00388
00389 irplib_irstd * irplib_std_get_type(
00390 double ra,
00391 double dec)
00392 {
00393 char ** catalog_names ;
00394 const irplib_irstd * refstar ;
00395 int nfound ;
00396 const irplib_irstd ** refstars ;
00397 int i, j ;
00398
00399
00400 catalog_names = (char**)irplib_irlist_catalogs ;
00401 nfound = 0 ;
00402
00403
00404 for (i=0 ; catalog_names[i] ; i++) {
00405 irplib_std_setactive("none");
00406 irplib_std_setactive(catalog_names[i]) ;
00407 refstar = irplib_std_get_closest_star(ra, dec) ;
00408 if (refstar != NULL) nfound ++ ;
00409 }
00410 refstars = cpl_malloc(nfound*sizeof(irplib_irstd*)) ;
00411 j = 0 ;
00412 for (i=0 ; catalog_names[i] ; i++) {
00413 irplib_std_setactive("none");
00414 irplib_std_setactive(catalog_names[i]) ;
00415 refstar = irplib_std_get_closest_star(ra, dec) ;
00416 if (refstar != NULL) {
00417 refstars[j] = refstar ;
00418 j++ ;
00419 }
00420 }
00421
00422
00423 refstar = NULL ;
00424 for (i=0 ; i<nfound ; i++) {
00425 if (strcmp(refstars[i]->sptype, "--")) {
00426 refstar = refstars[i] ;
00427 break ;
00428 }
00429 }
00430
00431 cpl_free(refstars) ;
00432 return (irplib_irstd*)refstar ;
00433 }
00434
00435
00445
00446 irplib_irstd * irplib_std_get_mag(
00447 double ra,
00448 double dec,
00449 irplib_band band,
00450 double * mag)
00451 {
00452 char ** catalog_names ;
00453 const irplib_irstd * refstar ;
00454 int nfound ;
00455 const irplib_irstd ** refstars ;
00456 int i, j ;
00457
00458
00459 catalog_names = (char**)irplib_irlist_catalogs ;
00460 nfound = 0 ;
00461
00462
00463 for (i=0 ; catalog_names[i] ; i++) {
00464 irplib_std_setactive("none");
00465 irplib_std_setactive(catalog_names[i]) ;
00466 refstar = irplib_std_get_closest_star(ra, dec) ;
00467 if (refstar != NULL) nfound ++ ;
00468 }
00469 refstars = cpl_malloc(nfound*sizeof(irplib_irstd*)) ;
00470 j = 0 ;
00471 for (i=0 ; catalog_names[i] ; i++) {
00472 irplib_std_setactive("none");
00473 irplib_std_setactive(catalog_names[i]) ;
00474 refstar = irplib_std_get_closest_star(ra, dec) ;
00475 if (refstar != NULL) {
00476 refstars[j] = refstar ;
00477 j++ ;
00478 }
00479 }
00480
00481
00482 refstar = NULL ;
00483 switch (band) {
00484 case BAND_J:
00485 for (i=0 ; i<nfound ; i++) {
00486 if (refstars[i]->mag_J < 98.0) {
00487 refstar = refstars[i] ;
00488 *mag = (double)(refstar->mag_J) ;
00489 break ;
00490 }
00491 }
00492 break ;
00493 case BAND_H:
00494 for (i=0 ; i<nfound ; i++) {
00495 if (refstars[i]->mag_H < 98.0) {
00496 refstar = refstars[i] ;
00497 *mag = (double)(refstar->mag_H) ;
00498 break ;
00499 }
00500 }
00501 break ;
00502 case BAND_K:
00503 for (i=0 ; i<nfound ; i++) {
00504 if (refstars[i]->mag_K < 98.0) {
00505 refstar = refstars[i] ;
00506 *mag = (double)(refstar->mag_K) ;
00507 break ;
00508 }
00509 }
00510 break ;
00511 case BAND_KS:
00512 for (i=0 ; i<nfound ; i++) {
00513 if (refstars[i]->mag_Ks < 98.0) {
00514 refstar = refstars[i] ;
00515 *mag = (double)(refstar->mag_Ks) ;
00516 break ;
00517 }
00518 }
00519 break ;
00520 case BAND_L:
00521 for (i=0 ; i<nfound ; i++) {
00522 if (refstars[i]->mag_L < 98.0) {
00523 refstar = refstars[i] ;
00524 *mag = (double)(refstar->mag_L) ;
00525 break ;
00526 }
00527 }
00528 break ;
00529 case BAND_M:
00530 for (i=0 ; i<nfound ; i++) {
00531 if (refstars[i]->mag_M < 98.0) {
00532 refstar = refstars[i] ;
00533 *mag = (double)(refstar->mag_M) ;
00534 break ;
00535 }
00536 }
00537 break ;
00538 case BAND_LP:
00539 for (i=0 ; i<nfound ; i++) {
00540 if (refstars[i]->mag_Lp < 98.0) {
00541 refstar = refstars[i] ;
00542 *mag = (double)(refstar->mag_Lp) ;
00543 break ;
00544 }
00545 }
00546 break ;
00547 case BAND_MP:
00548 for (i=0 ; i<nfound ; i++) {
00549 if (refstars[i]->mag_Mp < 98.0) {
00550 refstar = refstars[i] ;
00551 *mag = (double)(refstar->mag_Mp) ;
00552 break ;
00553 }
00554 }
00555 break ;
00556 default:
00557 break ;
00558 }
00559
00560
00561
00562 cpl_free(refstars) ;
00563 return (irplib_irstd*)refstar ;
00564 }
00565
00566
00572
00573 const char * irplib_std_band_name(irplib_band band)
00574 {
00575 switch (band) {
00576 case BAND_J: return "J" ;
00577 case BAND_JS: return "Js" ;
00578 case BAND_JBLOCK: return "J+Block" ;
00579 case BAND_H: return "H" ;
00580 case BAND_K: return "K" ;
00581 case BAND_KS: return "Ks" ;
00582 case BAND_L: return "L" ;
00583 case BAND_M: return "M" ;
00584 case BAND_LP: return "Lp" ;
00585 case BAND_MP: return "Mp" ;
00586 case BAND_Z: return "Z" ;
00587 case BAND_SZ: return "SZ" ;
00588 case BAND_SH: return "SH" ;
00589 case BAND_SK: return "SK" ;
00590 case BAND_SL: return "SL" ;
00591 default: return "Unknown" ;
00592 }
00593 }
00594
00595
00605
00606 const char * irplib_std_catalog_name(int cat_id)
00607 {
00608 return irplib_irlist_catalogs[cat_id] ;
00609 }
00610
00611
00612
00619
00620 const char * irplib_std_get_name(const irplib_irstd * self)
00621 {
00622
00623 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, NULL);
00624
00625 return self->name;
00626 }
00627
00628
00635
00636 const char * irplib_std_get_type_spectral(const irplib_irstd * self)
00637 {
00638
00639 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, NULL);
00640
00641 return self->sptype;
00642 }
00643
00644
00645
00652
00653 const char * irplib_std_get_catalog(const irplib_irstd * self)
00654 {
00655
00656 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, NULL);
00657
00658 return irplib_irlist_catalogs[self->source];
00659 }
00660
00661
00664
00678
00679 static int irplib_std_setactive(const char * catalog)
00680 {
00681 int i ;
00682 int found ;
00683 int active ;
00684
00685
00686 if (catalog==NULL) {
00687 found=0 ;
00688 i=0 ;
00689 while (irplib_irstd_list[i].name != NULL) {
00690 if (irplib_irstd_list[i].select) found++ ;
00691 i++ ;
00692 }
00693 return found ;
00694 }
00695
00696
00697 if (!strcmp(catalog, "none")) {
00698 i=0 ;
00699 while (irplib_irstd_list[i].name!=NULL) {
00700 irplib_irstd_list[i].select=0 ;
00701 i++;
00702 }
00703 return 0 ;
00704 }
00705
00706
00707 if (!strcmp(catalog, "all")) {
00708 i=0 ;
00709 while (irplib_irstd_list[i].name!=NULL) {
00710 irplib_irstd_list[i].select=1 ;
00711 i++;
00712 }
00713 return i ;
00714 }
00715
00716
00717 i=0 ;
00718 active=0 ;
00719 found=0 ;
00720 while (irplib_irstd_list[i].name!=NULL) {
00721 if (!strcmp(catalog, irplib_irlist_catalogs[irplib_irstd_list[i].source])) {
00722 found=1 ;
00723 irplib_irstd_list[i].select=1 ;
00724 }
00725 if (irplib_irstd_list[i].select) {
00726 active++ ;
00727 }
00728 i++ ;
00729 }
00730 if (found<1) return -1 ;
00731 return active ;
00732 }
00733
00734 #define IRPLIB_STD_MAXRADIUS (2.0/60.0)
00735 #define IRPLIB_STD_SQMAXRADIUS (IRPLIB_STD_MAXRADIUS*IRPLIB_STD_MAXRADIUS)
00736
00749
00750 static const irplib_irstd * irplib_std_get_closest_star(
00751 double ra_d,
00752 double dec_d)
00753 {
00754 const irplib_irstd * starlist ;
00755 double ra, dec ;
00756 double cur_dist ;
00757 double min_dist ;
00758 int min_index ;
00759 int i ;
00760
00761 starlist = NULL ;
00762
00763
00764 i=0 ;
00765 while (irplib_irstd_list[i].select==0 && irplib_irstd_list[i].name!=NULL)
00766 i++ ;
00767 if (irplib_irstd_list[i].name==NULL) return NULL ;
00768 ra = irplib_irstd_list[i].ra ;
00769 dec = irplib_irstd_list[i].dec ;
00770 min_dist = (ra_d-ra)*(ra_d-ra)+(dec_d-dec)*(dec_d-dec) ;
00771 min_index = i ;
00772
00773
00774 i=0 ;
00775 while (irplib_irstd_list[i].name!=NULL) {
00776 if (irplib_irstd_list[i].select) {
00777 ra = irplib_irstd_list[i].ra ;
00778 dec = irplib_irstd_list[i].dec ;
00779 cur_dist = (ra_d-ra)*(ra_d-ra)+(dec_d-dec)*(dec_d-dec) ;
00780 if (cur_dist<min_dist) {
00781 min_dist = cur_dist ;
00782 min_index = i ;
00783 }
00784 }
00785 i++ ;
00786 }
00787
00788 if (min_dist>IRPLIB_STD_SQMAXRADIUS) {
00789
00790 starlist = NULL ;
00791 } else {
00792 starlist = &(irplib_irstd_list[min_index]);
00793 }
00794 return starlist ;
00795 }
00796
00797