home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / fblnk224.zip / noninteractive.c < prev    next >
C/C++ Source or Header  |  1999-01-28  |  12KB  |  415 lines
  1. /*  noninteractive.c  */
  2. /*  part of the fitsblink program  */
  3. /*  Jure Skvarc, July 1998  */
  4. #include <stdio.h>
  5. #include <stdlib.h>
  6. #include <math.h>
  7. #include <string.h>
  8. #include <sys/types.h> 
  9. #include <sys/time.h>
  10. #include <unistd.h>
  11. #include <getopt.h>
  12. #include <time.h>
  13.  
  14. #include "functs.h"
  15. #include "consts.h"
  16.  
  17. extern STATE state;
  18. extern BLINK_FRAME *frame;
  19. extern CATALOG catalog[NUMBER_OF_CATALOGS];
  20. extern OBSERVATORY observatory;
  21. extern MAIL mail;
  22. extern FITS_IMAGE count_options;
  23. extern MATCH_OPTIONS match_options;
  24. extern TRANSROT trglob;
  25. extern char *optarg;
  26. extern int optopt, optind, opterr;
  27. /*  Measure star positions and output them to a file in non-interactive mode  */
  28.  
  29. void
  30. starcount(int argc, char *argv[])
  31.  
  32. {
  33.   int i;
  34.   int opt;
  35.  
  36.   while ((opt = getopt(argc, argv, "l:s:m:i:c:b:n:o:agv")) != -1) {
  37.     switch(opt) {
  38.     case 'l':
  39.       /*  Name of the log file  */
  40.       state.logfile = strdup(optarg);
  41.       if ((state.logfp = fopen(state.logfile,"a") )== NULL) {
  42.     fprintf(stderr, "Can not open log file!\n");
  43.     exit(1);
  44.       }
  45.       break;
  46.     case 'm':
  47.       count_options.maxsize = strtol(optarg, NULL, 10);
  48.       if (count_options.maxsize < 10) count_options.maxsize = 10;
  49.       break;
  50.     case 'b':
  51.       count_options.border = strtol(optarg, NULL, 10);
  52.       break;
  53.     case 'c':
  54.       count_options.minstar = strtol(optarg, NULL, 10);
  55.       break;
  56.     case 's':
  57.       /*  Sensitivity expressed in number of standard deviations above the 
  58.           background  */
  59.       count_options.sg_num = strtod(optarg, NULL);
  60.       if (count_options.sg_num < 0) count_options.sg_num = 2.0;
  61.       break;
  62.     case 'i':
  63.       /*  The minimum acceptable intensity  */
  64.       count_options.minbright = strtod(optarg, NULL);
  65.       break;
  66.     case 'n':
  67.       /*  The radius of the inner circle for the aperture astrometry  */
  68.       count_options.inner = strtod(optarg, NULL);
  69.       break;
  70.     case 'o':
  71.       /*  The radius of the outer circle for the aperture astrometry  */
  72.       count_options.outer = strtod(optarg, NULL);
  73.       break;
  74.     case 'a':
  75.       /*  Aperture astrometry  */
  76.       count_options.aperture = 1;
  77.       break;
  78.     case 'g':
  79.       /*  "Above sigma astrometry" */
  80.       count_options.aperture = 0;
  81.       break;
  82.     case 'v':
  83.       /*  Verbose mode */
  84.       state.verbose = 1;
  85.       break;
  86.     case '?':
  87.       switch(optopt) {
  88.       default:
  89.     fprintf(stderr, "Unknown option -%c ignored!\n", optopt);
  90.     break;
  91.       }
  92.       break;
  93.     }
  94.   }
  95.   for (i = optind; i < argc; i++) {
  96.     frame[0].name = strdup(argv[i]);
  97.     load_frame(frame, 0);
  98.     detect_stars(frame, 0);
  99.     star_list(frame, 0);
  100.     remove_frame(frame, 0);
  101.   }
  102. }
  103.  
  104.  
  105.  
  106. /*  Update WCS values in the FITS file  */
  107. int
  108. update_wcs(BLINK_FRAME *frame, int n, TRANSROT tr)
  109.  
  110. {  
  111.   fitsfile *fptr;
  112.   int err = 0;
  113.   double xinc, yinc;
  114.   double ac, dc, xc, yc, nw, nh, fi;
  115.  
  116.   /*  Open file for writing  */
  117.   fits_open_file(&fptr, frame[n].name, READWRITE, &err);
  118.   if (err != 0) {
  119.     if (!state.interactive) {
  120.       fprintf(stderr, "Can't open %s for writing!\n", frame[n].name);
  121.       return 1;
  122.     }
  123.   }
  124.   /*  Calculate WCS values  */
  125.   xinc = frame[n].map.width / frame[n].map.w;
  126.   yinc = frame[n].map.height / frame[n].map.h;
  127.   xc = 0.5 * frame[n].map.w;
  128.   yc = 0.5 * frame[n].map.h;
  129.   cal_new_wcs(tr, frame[n].map.ra, frame[n].map.dec, frame[n].map.fi, 
  130.           xc, yc, xinc, yinc, 
  131.           &ac, &dc, &nw, &nh, &fi);
  132.   fits_update_key(fptr, TSTRING, "CTYPE1", "RA---TAN", "Inserted by fitsblink", &err);
  133.   if (err != 0) {
  134.     fprintf(stderr, "Disaster: could not update keyword CTYPE1 for file %s\n", frame[n].name);
  135.     return 1;
  136.   }
  137.   fits_update_key(fptr, TDOUBLE, "CRPIX1", &xc, "Inserted by fitsblink", &err);
  138.   if (err != 0) {
  139.     fprintf(stderr, "Disaster: could not update keyword CRPIX1 for file %s\n", frame[n].name);
  140.     return 1;
  141.   }
  142.   fits_update_key(fptr, TDOUBLE, "CRVAL1", &ac, "Inserted by fitsblink", &err);
  143.   if (err != 0) {
  144.     fprintf(stderr, "Disaster: could not update keyword CRVAL1 for file %s\n", frame[n].name);
  145.     return 1;
  146.   }
  147.   fits_update_key(fptr, TDOUBLE, "CDELT1", &nw, "Inserted by fitsblink", &err);
  148.   if (err != 0) {
  149.     fprintf(stderr, "Disaster: could not update keyword CDELT1 for file %s\n", frame[n].name);
  150.     return 1;
  151.   }
  152.   fits_update_key(fptr, TDOUBLE, "CROTA1", &fi, "Inserted by fitsblink", &err);
  153.   if (err != 0) {
  154.     fprintf(stderr, "Disaster: could not update keyword CROTA1 for file %s\n", frame[n].name);
  155.     return 1;
  156.   }
  157.   fits_update_key(fptr, TSTRING, "CTYPE2", "DEC--TAN", "Inserted by fitsblink", &err);
  158.   if (err != 0) {
  159.     fprintf(stderr, "Disaster: could not update keyword CTYPE2 for file %s\n", frame[n].name);
  160.     return 1;
  161.   }
  162.   fits_update_key(fptr, TDOUBLE, "CRPIX2", &yc, "Inserted by fitsblink", &err);
  163.   if (err != 0) {
  164.     fprintf(stderr, "Disaster: could not update keyword CRPIX2 for file %s\n", frame[n].name);
  165.     return 1;
  166.   }
  167.   fits_update_key(fptr, TDOUBLE, "CRVAL2", &dc, "Inserted by fitsblink", &err);
  168.   if (err != 0) {
  169.     fprintf(stderr, "Disaster: could not update keyword CRVAL2 for file %s\n", frame[n].name);
  170.     return 1;
  171.   }
  172.   fits_update_key(fptr, TDOUBLE, "CDELT2", &nh, "Inserted by fitsblink", &err);
  173.   if (err != 0) {
  174.     fprintf(stderr, "Disaster: could not update keyword CDELT2 for file %s\n", frame[n].name);
  175.     return 1;
  176.   }
  177.   fits_update_key(fptr, TDOUBLE, "CROTA2", &fi, "Inserted by fitsblink", &err);
  178.   if (err != 0) {
  179.     fprintf(stderr, "Disaster: could not update keyword CROTA2 for file %s\n", frame[n].name);
  180.     return 1;
  181.   }
  182.   fits_close_file(fptr, &err);
  183.   if (err != 0) {
  184.     fprintf(stderr, "Disaster: could not close file %s\n", frame[n].name);
  185.     return 1;
  186.   }
  187.   return 0;
  188. }
  189.  
  190. /*  Read a star list, field center coordinates and pixel size (or fits header), */
  191. /*  compare it to catalog and output celestial coordinates along with any  */
  192. /*  cross ID's                                                             */
  193. void
  194. make_catalog(int argc, char *argv[])
  195.  
  196. {
  197.   int opt;
  198.   int wcs = 0;
  199.   char *starlist = NULL;
  200.   double ra, dec, p1, p2;
  201.   FILE *fp;
  202.   char a[256], fn[256];
  203.   char *cat = NULL, *path;
  204.   int nstar, n, i;
  205.   int xw, yw;
  206.   int writewcs = 0;
  207.  
  208.   match_options.firstbright = FIRSTBRIGHT;
  209.   match_options.nstar = CONSTSTAR;
  210.   match_options.min = CONSTMIN;
  211.   match_options.poserrmax = POSERRMAX;
  212.   match_options.magerrmax = MAGERRMAX;
  213.   match_options.maxres = MAXRES;
  214.   
  215.   while ((opt = getopt (argc, argv, "b:l:n:d:a:w:h:m:e:i:r:f:c:x:y:vW")) != -1) {
  216.     switch(opt) {
  217.     case 'l':
  218.       /*  Name of the log file  */
  219.       state.logfile = strdup(optarg);
  220.       if ((state.logfp = fopen(state.logfile,"a") )== NULL) {
  221.     fprintf(stderr, "Can not open log file!\n");
  222.     exit(1);
  223.       }
  224.       break;
  225.     case 'b':
  226.       match_options.firstbright = strtol(optarg, NULL, 10);
  227.       if (match_options.firstbright < MINFIRSTBRIGHT) match_options.firstbright = MINFIRSTBRIGHT;
  228.       break;
  229.     case 'a':
  230.       /* Get right ascension  */
  231.       ra = strtod(optarg, NULL);
  232.       wcs++;
  233.       break;
  234.     case 'd':
  235.       /* Get declination */
  236.       dec = strtod(optarg, NULL);
  237.       wcs++;
  238.       break;
  239.     case 'w':
  240.       /* Get pixel width */
  241.       p1 = strtod(optarg, NULL);
  242.       wcs++;
  243.       break;
  244.     case 'h':
  245.       /* Get pixel height */
  246.       p2 = strtod(optarg, NULL);
  247.       wcs++;
  248.       break;
  249.     case 'x':
  250.       /*  Number of pixels in x direction  */
  251.       xw = strtol(optarg, NULL, 10);
  252.       wcs++;
  253.       break;
  254.     case 'y':
  255.       /*  Number of pixels in y direction  */
  256.       yw = strtol(optarg, NULL, 10);
  257.       wcs++;
  258.       break;
  259.     case 'f':
  260.       /*  Name of the star list file  */
  261.       starlist = strdup(optarg);
  262.       if ((fp = fopen(starlist, "r") )== NULL) {
  263.     fprintf(stderr, "Can not open starlist file: %s!\n", starlist);
  264.     exit(1);
  265.       }
  266.     case 'n':
  267.       match_options.nstar = strtol(optarg, NULL, 10);
  268.       break;
  269.     case 'm':
  270.       match_options.min = strtod(optarg, NULL);
  271.       break;
  272.     case 'e':
  273.       match_options.poserrmax = strtod(optarg, NULL);
  274.       break;
  275.     case 'i':
  276.       match_options.magerrmax = strtod(optarg, NULL);
  277.       break;
  278.     case 'r':
  279.       match_options.maxres = strtod(optarg, NULL);
  280.       break;
  281.     case 'c':
  282.       /*  Extract catalog name and path  */
  283.       cat = optarg;
  284.       path = strtok(cat, ":");
  285.       path = strtok(NULL, ":");
  286.       break;
  287.     case 'v':
  288.       state.verbose = 1;
  289.       break;
  290.     case 'W':
  291.       writewcs = 1;
  292.       break;
  293.     case '?':
  294.       switch(optopt) {
  295.       default:
  296.     fprintf(stderr, "Unknown option -%c ignored!\n", optopt);
  297.     break;
  298.       }
  299.       break;
  300.     }
  301.   }
  302.  
  303.   if (starlist == 0) {
  304.     fprintf(stderr, "No file name!\n");
  305.     exit(1);
  306.   }
  307.   /*  In the first file line we look for the image file name  */
  308.   fscanf(fp, "Input file name: %s", fn);
  309.   frame[0].name = strdup(fn);
  310.   /*  Load the image  */
  311.   if (load_frame(frame, 0)) {
  312.     fprintf(stderr, "Can't load %s\n", frame[0].name);
  313.     exit(1);
  314.   }
  315.   if (wcs != 6) {
  316.     if (!frame[0].fitsfile.wcs) {
  317.       fprintf(stderr, "No image coordinates!\n");
  318.       exit(1);
  319.     }
  320.   }
  321.   else {
  322.     if (!frame[0].fitsfile.wcs) {
  323.       frame[0].fitsfile.nax[0] = xw;
  324.       frame[0].fitsfile.nax[1] = yw;
  325.     }
  326.     /*  Write or override image coordinates and pixel dimensions  */
  327.     frame[0].ra = frame[0].fitsfile.crval1 = ra;
  328.     frame[0].dec = frame[0].fitsfile.crval2 = dec;
  329.     frame[0].fitsfile.crpix1 = 0.5 * frame[0].fitsfile.nax[0];
  330.     frame[0].fitsfile.crpix2 = 0.5 * frame[0].fitsfile.nax[1];
  331.     frame[0].xsize = p1; 
  332.     frame[0].ysize = p2; 
  333.     frame[0].fitsfile.cdelt1 = p1 / 3600.0;
  334.     frame[0].fitsfile.cdelt2 = p2 / 3600.0;
  335.     frame[0].fitsfile.crot = 0;
  336.     strcpy(frame[0].fitsfile.ctype, "-TAN");
  337.     frame[0].fitsfile.wcs = 1;
  338.   }
  339.   /*  reserve some space for the star list  */
  340.   nstar = NSTAR;
  341.   frame[0].imagelist.s = (STAR *) myalloc(sizeof(STAR) * nstar, "make_catalog", 
  342.                       "frame[0].star.s");
  343.   /*  Read a star list  */
  344.   /*====================*/
  345.   i = 0;
  346.   while (!feof(fp)) {
  347.     fgets(a, 256, fp);    /*  Read a complete line   */
  348.     n = sscanf(a, "%d %lf %lf %f", &frame[0].imagelist.s[i].n,
  349.        &frame[0].imagelist.s[i].x, &frame[0].imagelist.s[i].y,
  350.        &frame[0].imagelist.s[i].mag);
  351.     if (n == 4 && !feof(fp)) {
  352.       i++;
  353.       if (i == nstar) {
  354.     nstar += 250;
  355.     frame[0].imagelist.s = (STAR *) 
  356.       realloc(frame[0].imagelist.s, sizeof(STAR) * nstar);
  357.     if (frame[0].imagelist.s == NULL) {
  358.       fprintf(stderr, "Too many stars!");
  359.       exit(1);
  360.     }
  361.       }
  362.     }
  363.   }
  364.   frame[0].imagelist.n = i;
  365.   fclose(fp);
  366.   if (cat == NULL) {
  367.     fprintf(stderr, "Missing catalog name!\n");
  368.     exit(1);
  369.   }
  370.   /*  Check which catalog we have  */
  371.   /*===============================*/
  372.   if (strcmp(cat, "usno") == 0) {
  373.     catalog[CAT_USNO].use = 1;
  374.     catalog[CAT_USNO].path = strdup(path);
  375.     catalog[CAT_USNO].lpath = strlen(path);
  376.  
  377.   } 
  378.   else if (strcmp(cat, "gscn") == 0) {
  379.     catalog[CAT_GUIDE].use = 1;
  380.     catalog[CAT_GUIDE].path = strdup(path);
  381.     catalog[CAT_GUIDE].lpath = strlen(path);
  382.   } 
  383.   else if (strcmp(cat, "gscs") == 0) {
  384.     catalog[CAT_GUIDE_SOUTH].use = 1;
  385.     catalog[CAT_GUIDE_SOUTH].path = strdup(path);
  386.     catalog[CAT_GUIDE_SOUTH].lpath = strlen(path);
  387.   } 
  388.   /*  else if (strcmp(cat, "ppm") == 0) {
  389.     catalog[CAT_PPM].use = 1;
  390.     catalog[CAT_PPM].path = strdup(path);
  391.     catalog[CAT_PPM].lpath = strlen(path);
  392.     } */
  393.   else {
  394.     fprintf(stderr, "Unknown catalog\n");
  395.     fprintf(stderr, "Try one of these:\n");
  396.     fprintf(stderr, "  usno -- USNO SA 1.0\n");
  397.     fprintf(stderr, "  gscn -- GSC 1.1 from 90 to -7.5 declination\n");
  398.     fprintf(stderr, "  gscs -- GSC 1.1 from -7.5 to -90 declination\n");
  399.     /*    fprintf(stderr, "  ppm -- PPM\n"); */
  400.     exit(1);
  401.   }
  402.   frame[0].match = match_options;
  403.   /*  Match a list with a star catalog  */
  404.   if (match_stars(frame, 0)) {
  405.     fprintf(stderr, "Matching failed!\n");
  406.     exit(1);
  407.   }
  408.   /*  Output star list  */
  409.   star_list(frame, 0);
  410.   if (writewcs) {
  411.     update_wcs(frame, 0, trglob);
  412.   }
  413.   remove_frame(frame, 0);
  414. }
  415.