/** $Header: /green/st2/src/iros/RCS/cerr.c,v 1.7 1998/09/23 15:37:41 jaaps Exp $ ********************************** - - **************************************** * * _/_/_/_/_/_/ _/_/_/ _/_/ _/_/ SRON Utrecht * _/_/_/_/ _/_/ _/_/_/_/_/ _/_/_/ _/_/_/ * _/_/_/ _/_/ _/_/ _/_/_/_/_/ _/ _/ _/_/_/_/ _/_/_/ * _/_/_/_/_/ _/_/_/_/_/_/_/ _/_/_/ _/ _/ _/ _/_/ * _/_/_/_/_/_/_/_/_/_/ _/_/_/_/_/ _/_/ _/ _/ _/ _/_/_/ _/ * _/_/ _/_/_/_/_/_/ _/_/_/_/_/ _/_/_/ _/_/ _/_/ _/ _/_/ * _/_/_/_/_/_/ _/_/ _/_/_/_/ _/_/ _/ _/ _/ _/_/_/ * * - - - - - - - - - - - - - - - - - - - - - - - * Use : cerr -p .resfile * * Parameters : resfile: an iros outputfile, "*.res". * * - - - - - - - - - - - - - - - - - - - - - - - * Dependencies: * * - - - - - - - - - - - - - - - - - - - - - - - * Description : report some statistics on the deviations in the source * positions as listed in the catalog and those found by * iros. * * - - - - - - - - - - - - - - - - - - - - - - - * $Log: cerr.c,v $ * Revision 1.7 1998/09/23 15:37:41 jaaps * check if the catalog file can be opened and read properly. * before this, open/read failures produced a message stating * that no sources could be identified. * * Revision 1.6 1998/09/22 12:00:11 jaaps * altered d_angle computation, problems with slalib installations. * * Revision 1.5 1998/06/18 07:42:08 jaaps * non essential, cosmetic adjustments (no reason found in source). * * Revision 1.4 1998/05/13 10:28:32 jaaps * non essential, cosmetic adjustments (no reason found in source). * * Revision 1.3 1997/11/12 14:39:09 jaaps * src_readcatalog parameter list adapted, firts arg. now existing cataloglist or NULL. * fixed some leftover stuff and loose ends. * suppress use of source-entries in those energybands where it obviously * was not detected, obvious here has to be 0 flux and 0 fluxeror. * * Revision 1.2 1997/09/22 13:37:24 jaaps * non essential, cosmetic adjustments (no reason found in source). * * Revision 1.1 1997/09/08 10:08:59 jaaps * Initial revision * * ********************************** - - ****************************************/ #include #include #include #include #include #include "wfcglobals.h" #include "instrument.h" #include "expose.h" #include "iros.h" #include "sources.h" #include "util.h" #define C_CALLOC(s,a,b) secure_alloc (#s, a, b) /** *** #define C_CALLOC(s,a,b) calloc (a, b) **/ typedef struct { float epoch; double ra, dec; } coordinate; typedef struct { cardinal npln; TiLwHh *bands; cardinal *nn; float *xsum, *ysum, *xsumsq, *ysumsq, *ts, *tsq; } STATISTICSREC, *STATISTICSPTR; static char *RCSID= "$Header: /green/st2/src/iros/RCS/cerr.c,v 1.7 1998/09/23 15:37:41 jaaps Exp $"; int gbl_debug = 0; char *gbl_creator= "$State: Exp $"; #ifndef PACKET_REVN2 #define PACKET_REVN2 9999 #endif int gbl_revision= PACKET_REVN2; float CCFOR(sla_dsep) (double *, double *, double *, double *); float d_angle (char *, char *, coordinate *, coordinate *); typedef struct { char *resfile; int mode; float fl, fclln, fov; rpair align; } INPPARREC; static KeyRecord par_key[]= { {"RESFILE" , FTP(pchar, resfile, INPPARREC, PFC_IFILE), 0, "input, iros (.res) resultfile"}, {"FLUXLIM" , FTP(real , fl , INPPARREC, PFC_DFLT ), 0, "lower limit on sourceflux."}, {"XY" , FTP(int , mode , INPPARREC, PFC_DFLT ), 0, "set 1 to report pixel diferences"}, {"PALIGN" , FTP(rpair, align , INPPARREC, PFC_DFLT ), 0, ".res pointing update, x,y"}, {"FLDOFVW" , FTP(rpair, fov , INPPARREC, PFC_DFLT ), 0, "field of view (degrees)"}, {"FOCALLNG", FTP(real , fclln , INPPARREC, PFC_DFLT ), 0, "overrides instrument focallength"} }; #define npar_key (sizeof (par_key) / sizeof (par_key[0])) PARSW2KEYREC swmap[]= { "a" , "PALIGN", "F" , "FOCALLNG", "l" , "FLUXLIM", "V" , "FLDOFVW", "x" , "XY", NULL , NULL }; char *reqArg[]= {"RESFILE", NULL}; cardinal NreqArg = 1; cardinal get_cat_entry (char *name, WFC_CATALOGPTR ct) { cardinal i; for (i= 0; i!= ct->n; i++) if ((strcmp (name, wft_btAget (ct->e, i, CR_CNAME)) == 0)) return i; return ct->n; } void d_header0 (void) { wLog (DBG_CHAT, " catalog-position iros-position diff (arcmin)\n"); } void d_trailer0 (STATISTICSPTR sts, float fl) { char mess[256]; float mn= 0.0, sg= 0.0, px= 0.0; if (sts->nn[0] > 0) mn= sts->xsum[0]/(float) sts->nn[0]; if (sts->nn[0] > 1) sg= (sts->xsumsq[0]-sts->xsum[0]*mn)/(float) (sts->nn[0]-1); px= fl*tan (WFC_DEG2RAD(mn)); sprintf (mess, "\n mean difference %.3f (%.3f) arcmin (== %.3f pix.)\n", 60.0*mn, 60.0*sg, px); wLog (DBG_CHAT, mess); } void d_header1 (void) { wLog (DBG_CHAT, " catalog position iros position\n"); wLog (DBG_CHAT, " name x y chnls x y dx dy length\n"); wLog (DBG_CHAT, "\n"); } void d_trailer1 (STATISTICSPTR sts) { char mess[256]; cardinal i, j; float t; wLog (DBG_CHAT, "\n"); wLog (DBG_CHAT, " E mean (pixels) length (pixels) angle (deg)\n"); wLog (DBG_CHAT, "chnl dx (s.d) dy (s.d) ||.|| (s.d)\n"); for (j=0; j!= sts->npln; j++) { if ((t= (float) sts->nn[j]) > 0.0) sprintf (mess, "%2.2d-%2.2d %8.3f (%5.2f) %8.3f (%5.2f) %8.3f (%5.2f) %6.1f\n", sts->bands[j].L, sts->bands[j].H, (sts->nn[j] > 0) ? sts->xsum[j]/t : 0.0, (sts->nn[j] > 1) ? (sts->xsumsq[j]-sts->xsum[j]*sts->xsum[j]/t)/(t-1.0) : 0.0, (sts->nn[j] > 0) ? sts->ysum[j]/t : 0.0, (sts->nn[j] > 1) ? (sts->ysumsq[j]-sts->ysum[j]*sts->ysum[j]/t)/(t-1.0) : 0.0, (sts->nn[j] > 0) ? sts->ts[j]/t : 0.0, (sts->nn[j] > 1) ? (sts->tsq[j]-sts->ts[j]*sts->ts[j]/t)/(t-1.0) : 0.0, WFC_RAD2DEG(atan2 (sts->ysum[j], sts->xsum[j]))); else sprintf (mess, "%2.2d-%2.2d ------- no data -------\n", sts->bands[j].L, sts->bands[j].H); wLog (DBG_CHAT, mess); } } float d_angle (char *name, char *tname, coordinate *cfrom, coordinate *ctows) { double d1, d0, a1, a0, na, phi, eta, N, sd0, sd1, cd0, cd1, dra, sdra, cdra; float rv; char t0[128], t1[128], mess[256]; a0 = WFC_DEG2RAD(cfrom->ra); d0 = WFC_DEG2RAD(cfrom->dec); a1 = WFC_DEG2RAD(ctows->ra); d1 = WFC_DEG2RAD(ctows->dec); /** ** used slaDs2tp form C-slalib here, link problems reported (SDC-rome) **/ sd0 = sin (d0); sd1 = sin (d1); cd0 = cos (d0); cd1 = cos (d1); dra = a0 - a1; sdra= sin (dra); cdra= cos (dra); N = sd0*sd1 + cd0*cd1*cdra; phi = cd0*dra/N; eta = (sd0*cd1 - cd0*sd1*cdra)/N; na= atan (sqrt (tan(phi)*tan(phi)+tan(eta)*tan(eta))); rv= (float) WFC_RAD2DEG(na); sprintf (mess, "%-10.10s (%-10.10s) %s %s %7.3f arcmin\n", (name == NULL || *name == '\0') ? "" : name, (tname == NULL || *tname == '\0') ? "" : tname, frmt_rd (cfrom->ra, cfrom->dec, t0), frmt_rd (ctows->ra, ctows->dec, t1), 60.0*rv); wLog (DBG_CHAT, mess); return rv; } STATISTICSPTR cerr_createstat (cardinal n) { STATISTICSPTR s= (STATISTICSPTR) C_CALLOC (cerr_createstat, 1, sizeof (STATISTICSREC)); if ((s->npln= n) > 0) { s->nn = (cardinal *) C_CALLOC (cerr_createstat, n, sizeof (cardinal)); s->bands = (TiLwHh *) C_CALLOC (cerr_createstat, n, sizeof (TiLwHh)); s->xsum = (float *) C_CALLOC (cerr_createstat, n, sizeof (float)); s->ysum = (float *) C_CALLOC (cerr_createstat, n, sizeof (float)); s->xsumsq= (float *) C_CALLOC (cerr_createstat, n, sizeof (float)); s->ysumsq= (float *) C_CALLOC (cerr_createstat, n, sizeof (float)); s->ts = (float *) C_CALLOC (cerr_createstat, n, sizeof (float)); s->tsq = (float *) C_CALLOC (cerr_createstat, n, sizeof (float)); } return s; } STATISTICSPTR cerr_radecdev (WFC_IROSPTR irs, STATISTICSPTR sts, float flim) { coordinate cc, ic; float rv; cardinal i, ck, maxplns; char mess[256], *name; if (sts == NULL) sts= cerr_createstat (1); for (i= 0; i!= irs->ctl->n; i++) { if (SOURCE_IS_BG(irs->ctl->e, i)) continue; name= wft_btAget (irs->ctl->e, i, CR_CNAME); if ((ck= get_cat_entry (name, irs->rsy->exp->ctl)) == irs->rsy->exp->ctl->n) { sprintf (mess, "cerr: cannot find source %s in %s\n", (name == NULL) ? "" : name, (irs->rsy->exp->ctl->name == NULL) ? "" : irs->rsy->exp->ctl->name); wLog (DBG_WARNING, mess); } else if (wft_btFget (irs->ctl->e, i, CR_FLUX, 0) >= flim) { ic.ra = wft_btDget (irs->ctl->e, i, CR_RA, 0); ic.dec= wft_btDget (irs->ctl->e, i, CR_DEC, 0); cc.ra = wft_btDget (irs->rsy->exp->ctl->e, ck, CR_RA, 0); cc.dec= wft_btDget (irs->rsy->exp->ctl->e, ck, CR_DEC, 0); rv= d_angle (wft_btAget (irs->rsy->exp->ctl->e, ck, CR_CNAME), wft_btAget (irs->rsy->exp->ctl->e, ck, CR_TNAME), &cc, &ic); sts->nn[0]++; sts->xsum[0] += rv; sts->xsumsq[0]+= rv*rv; } } return sts; } STATISTICSPTR cerr_xydev (WFC_IROSPTR irs, STATISTICSPTR sts, float fcln, float flim) { cardinal i, j, ck, maxplns, lmess; char *name, mess[256]; coordinate cc; float cx, cy, dx, dy, t; set_trafomatrix_ceim (&irs->rsy->exp->ra, &irs->rsy->exp->dec, &irs->rsy->exp->na); if (sts == NULL) sts= cerr_createstat (irs->rsy->exp->nplanes); maxplns= (sts->npln > irs->rsy->exp->nplanes) ? irs->rsy->exp->nplanes : sts->npln; for (i= 0; i!= irs->ctl->n; i++) { if (SOURCE_IS_BG(irs->ctl->e, i)) continue; name= wft_btAget (irs->ctl->e, i, CR_CNAME); if ((ck= get_cat_entry (name, irs->rsy->exp->ctl)) == irs->rsy->exp->ctl->n) { sprintf (mess, "cerr: cannot find source %s in %s\n", (name == NULL) ? "" : name, (irs->rsy->exp->ctl->name == NULL) ? "" : irs->rsy->exp->ctl->name); wLog (DBG_WARNING, mess); } else { cc.ra = wft_btDget (irs->rsy->exp->ctl->e, ck, CR_RA, 0); cc.dec= wft_btDget (irs->rsy->exp->ctl->e, ck, CR_DEC, 0); trafo_ceim (&cc.ra, &cc.dec, &cx, &cy, 1000.0); cx*= fcln; cy*= fcln; sprintf (mess, " %-10.10s %8.3f,%8.3f", wft_btAget (irs->ctl->e, i, CR_CNAME), cx, cy); lmess= strlen (mess); for (j=0; j!= maxplns; j++) { *mess= ' '; if (wft_btFget (irs->ctl->e, i, CR_FLUX, j) >= flim) { dx= cx-wft_btFget (irs->ctl->e, i, CR_X, j); dy= cy-wft_btFget (irs->ctl->e, i, CR_Y, j); sts->xsum[j] += dx; sts->ysum[j] += dy; sts->xsumsq[j]+= dx*dx; sts->ysumsq[j]+= dy*dy; t = sqrt ((double) (dx*dx+dy*dy)); sts->ts[j]+= t; sts->tsq[j]+= t*t; if (sts->nn[j] == 0) sts->bands[j]= irs->rsy->exp->bands[j]; if (j != 0) memset ((void *) mess, ' ', sizeof (mess)); sprintf (mess+lmess, " %2.2d-%2.2d %8.3f,%8.3f %8.3f,%8.3f %8.3f", sts->bands[j].L, sts->bands[j].H, wft_btFget (irs->ctl->e, i, CR_X, j), wft_btFget (irs->ctl->e, i, CR_Y, j), dx, dy, t); if (wft_btFget (irs->ctl->e, i, CR_DFLUX, j) == 0.0) { *mess= '*'; sprintf (mess+strlen (mess), " * 0 dflux!"); } wLog (DBG_CHAT, mess); sts->nn[j]++; } wLog (DBG_CHAT, "\n"); } } } return sts; } WFC_IROSPTR cerr_open (char *name) { WFC_IROSPTR irs; char mess[256]; irs= cor_open (name, TRUE); if (irs->ctl == NULL || irs->ctl->n == 0) { sprintf (mess, "cerr_open: no sourcelist in \"%s\"\n", name); wLog (DBG_WARNING, mess); return NULL; } sprintf (mess, "> \"%s\"\n", name); wLog (DBG_CHAT, mess); return irs; } void cerr_adjust (WFC_IROSPTR irs, float fl, rpair a) { char t0[256], mess[256]; float x, y; if (a.x != 0.0 || a.y != 0.0) { float x,y; sprintf (mess, "cerr: pointing update: %.4f, %.4f\n", a.x, a.y); wLog (DBG_CHAT, mess); sprintf (mess, " was: %s\n", frmt_rda (irs->rsy->exp->ra, irs->rsy->exp->dec, irs->rsy->exp->na, t0)); wLog (DBG_CHAT, mess); set_trafomatrix_imce (&irs->rsy->exp->ra, &irs->rsy->exp->dec, &irs->rsy->exp->na); x= a.x/fl; y= a.y/fl; trafo_imce (&x, &y, &irs->rsy->exp->ra, &irs->rsy->exp->dec); sprintf (mess, " becomes: %s\n\n", frmt_rda (irs->rsy->exp->ra, irs->rsy->exp->dec, irs->rsy->exp->na, t0)); wLog (DBG_CHAT, mess); } } main (int argc, char **argv) { INPPARREC par; FITS_HDU *hdu; WFC_IROSPTR irs; WFC_INSTRUMENTPTR ins; WFC_CATALOGPTR ctl; cardinal i, nfiles; char mess[256], **flist; STATISTICSPTR sts= NULL; float fov, fcln; memset ((void *) &par, 0, sizeof (INPPARREC)); hdu= new_hdu ((void *)&par, 1, 0, par_key, npar_key); wft_getparrev (argc, argv, NreqArg, reqArg, hdu, swmap, gbl_revision); free (hdu); hdu= NULL; fcln= fov= -1.0; if (par.mode == 0) d_header0 (); else d_header1 (); flist= wft_resolvelist (par.resfile, &nfiles); for (i= 0; i!= nfiles; i++) { if ((irs= cerr_open (flist[i])) != NULL) { if (fov <= 0.0) { ins= NULL; if (irs->rsy == NULL || irs->rsy->exp == NULL || irs->rsy->exp->insfile == NULL) { sprintf (mess, "cerr: no instrument ref. in res file \"%s\"\n", flist[i]); wLog (DBG_WARNING, mess); } else ins= instrument_get (irs->rsy->exp->insfile, NULL, 0); fcln= (ins != NULL) ? ins->focallength : par.fclln; fov = (ins != NULL) ? (0.5*sqrt((double) (ins->cor.sn.x*ins->cor.sn.x+ ins->cor.sn.y*ins->cor.sn.y)))/fcln : tan ((double) WFC_DEG2RAD (par.fov)); if (fcln <= 0.0 || fov <= 0.0) { sprintf (mess, "cerr: supply both mask-detector distance, param FOCALLNG (-F) and\n"); wLog (DBG_CHAT , mess); sprintf (mess, "---- the field of view, param FLDOFVW (-V), in degrees\n"); wLog (DBG_FATAL, mess); } ctl= src_readcatalog (NULL, irs->catfile, irs->rsy->exp->ra, irs->rsy->exp->dec, fov); if (ctl == NULL || ctl->n == 0) { sprintf (mess, "cerr: cannot read (or find any sources) in catalog %s\n", irs->catfile); wLog (DBG_WARNING, mess); sprintf (mess, " the pointing is %.3f,%.3f (ra,dec), search radius %.1f degr.\n", irs->rsy->exp->ra, irs->rsy->exp->dec, (double) WFC_RAD2DEG (atan((double) fov))); wLog (DBG_FATAL, mess); } } if (par.align.x != 0.0 || par.align.y != 0.0) cerr_adjust (irs, fcln, par.align); irs->rsy->exp->ctl= ctl; if (par.mode == 0) sts= cerr_radecdev (irs, sts, par.fl); else sts= cerr_xydev (irs, sts, fcln, par.fl); } } if (sts == NULL || sts->nn == NULL || sts->nn[0] == 0) { wLog (DBG_ERROR, "No sources could be identified (by name) in the catalog\n"); exit (0); } if (par.mode == 0) d_trailer0 (sts, fcln); else d_trailer1 (sts); }