class.c

Go to the documentation of this file.

/* class.c, write SMT data to a class data file */

#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <malloc.h>
#include <math.h>
#include <unistd.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>

#include "global.h"

#include "linuxlib_proto.h"


/* ---------------------------------------------------------------------- */
/*  ASTRO       Utility routine                                           */
/*      Returns Julian date of 1 day of Year                              */
/*  Arguments                                                             */
/*      NYR     I       Year                            Input             */
/*      JULDA   I       Date returned                   Output            */
/* ---------------------------------------------------------------------- */
int julda(int nyr)
{
  int nyrm1, ic, nadd, julday;
  int izero = 2025;

  nyrm1 = nyr - 1 ;

  ic = nyrm1 / 100;

  nadd = nyrm1 / 4 - ic + ic / 4;

  julday = 365. * (nyrm1 - izero) + nadd;

  return(julday);
}


/* ---------------------------------------------------------------------- */
/*  ASTRO       Utility routine                                           */
/*      Convert a date (id/im/iy) in number of days till the              */
/*      01/01/1984              (2025??)                                  */
/*  Arguments :                                                           */
/*      ID      I       Day number              Input                     */
/*      IM      I       Month number            Input                     */
/*      IY      I       Year number             Input                     */
/*      JOUR    I       Elapsed day count       Output                    */
/* ---------------------------------------------------------------------- */
void datj(int id,int im, int iy, int* jour)
{
  int ibiss, jour4;
  int ideb[13] = {0,31,59,90,120,151,181,212,243,273,304,334,365};

  if(im > 12 || im <= 0)
  {
    *jour = -99999;                                              /* Bad month */

    return;
  }

  jour4 = id + ideb[im-1];
  ibiss = julda(iy+1) - julda(iy) - 365;

  if(im >= 3)
    jour4 = jour4 + ibiss;

  if(im == 2)
  {
    if(id <= 0 || id > (ideb[im] - ideb[im-1] + ibiss)) 
    {
      *jour = -99999;                                     /* Bad February day */

      return;
    }
  }
  else
  {
    if(id <= 0 || id > (ideb[im] - ideb[im-1])) 
    {
      *jour = -99999;                                     /* Bad day of month */

      return;
    }
  }
  jour4 = jour4 + julda(iy);

  *jour = jour4;
}



                       /* convert pops date (yyyy.mmdd) into class day number */
int class_date(double pops_date)
{
  int day, month, year;
  double dd;
  static int jour;

  dd    = pops_date + 0.00005;                                    /* round up */
  year  = (int)dd;
  dd   -= (int)dd;
  dd   *= 100.0;
  month = (int)dd;
  dd   -= (int)dd;
  dd   *= 100.0;
  day   = (int)dd;


  datj(day, month, year, &jour);

  return(jour);
}


void add_class_section(long *info, int num, int type, int length)
{
  int i;

  if(num > 15)
  {
    st_fail("Scan header full, cannot add header section %d", type);
    return;
  }
                                                /* make space for new section */
  if(num == 1)
  {
    *(info + 2) = *info;
  }
  else
  {
    for(i= num - 2; i >= 0; i--)                            /* previous addrs */
    {
      *(info + 2 * num + i) = *(info + 2 * (num-1) + i);
    }
    for(i= num - 2; i >= 0; i--)                          /* previous lengths */
    {
      *(info + num + i) = *(info + (num-1) + i);
    }
  }
                                                          /* add section info */
  *(info + num - 1) = type;

  *(info + 2 * num - 1) = length;

  if(num > 1) 
  {
    *(info + 3 * num - 1) = *(info+ 3 * num - 2) + *(info + 2 * num - 2);
  }
}


              /* Convert (ra,dec) to galactic coordinats.  All are in degrees */
int eqtogal(double ra, double dec, double *glong, double *glat)
{
  double x, y, z, x0, z0, cosd;
  static double sin626=0, cos626;

  if(sin626 == 0)
  {
    sin626 = deg_sin(62.6);
    cos626 = deg_cos(62.6);
  }

  ra = ra - 192.25;
  cosd = deg_cos(dec);
  x0 = cosd * deg_cos(ra);
  x = cosd * deg_sin(ra);
  z0 = deg_sin(dec);
  y = -cos626 * x0 + sin626 * z0;
  z =  sin626 * x0 + cos626 * z0;

  *glat = (180./M_PI)*asin(z);
  *glong = (180./M_PI)*atan2(y,x) + 33;

  if(*glong < 0.0)
  {
    *glong += 360.0;
  }

  return(0);
}


int class_write(file, headp, datap, cname, verbose)
char *file, *cname;
struct HEADER *headp;
float *datap;
int verbose;
{
  int    fd, len_index, subscan;
  long   ent_off, scan_off, sect_off, data_off;
  char   sect_buf[REC_LEN];
  char  *cp, sourcename[256], srctmp[256], tmp[256], freqname[256], 
         tbuf[256], telescope[256];
  double glong_pos, glat_pos;
  struct scan_ent scan;
  struct file_index index;
  struct index_ent entry;
  float  scanFraction;
  union {
    struct general_descr general;
    struct position_descr position;
    struct spectro_descr spectro;
    struct fswitch_descr fswitch;
    struct calibration_descr calib;
  }sect;

  scanFraction = headp->scan - floor(headp->scan);

  if(verbose)
    st_report("class_write called");

  if(!strncmp(headp->obsmode, "LINEPCAL", 8))
  {
            /* don't record the normal cals, just the hot, cold and sky ones. */
    if(scanFraction < 0.205)
    {
      st_report("Obs %8.8s not processed for class, skip %7.2f.",
                headp->obsmode, headp->scan);

      return(0);
    }
  }

  fd = open(file, Oflag, Mode);

  if(fd == -1)
  {
    st_fail("Can't open class file: %s", file);

    return(0);
  }

  lseek(fd, 0L, SEEK_SET);
  len_index = read(fd, (char *) &index, sizeof (index));

                                                       /* is extension full? */
  if(((index.next_ent -1) != 0) && ((index.next_ent -1) % index.num_ent) == 0)
  {
    index.ext_addr[index.num_ext] = index.next_rec;
    index.next_rec += index.num_ent * ENT_LEN / REC_LEN;
    index.num_ext ++;
  }

        /* check if we handle this type of observation, then fill index entry */
  if(strncmp(headp->obsmode,"LINE", 4) == 0)
  {
    if(strncmp(&headp->obsmode[4], "OTF", 3) == 0)
    { 
      st_report("CLASS_WRITER: Obs type %8.8s not processed, skip scan %7.2f.",
        headp->obsmode, headp->scan);

      close(fd);

      return(0);
    }
    else
    {
      entry.obs_type = KIND_SPEC;
    }
  }
  else
  {
    st_fail("CLASS_WRITER: Unknown Obs type %8.8s, skip scan %7.2f.",
      headp->obsmode, headp->scan);

    close(fd);

    return(0);
  }

  entry.addr = index.next_rec;
  entry.obs = index.next_ent;
  index.next_ent ++;
  entry.version = 1;

  strncpy(srctmp, headp->object, 16);
  srctmp[16] = '\0';

  cp = strtok(srctmp, " ");
  if(cp)
  {
    sprintf(sourcename, "%-12.12s", cp);
  }
  else
  {
    strcpy(sourcename, "No Src Name ");
  }

  bcopy (sourcename, entry.source, 12);

  sprintf(freqname, "%-12s", headp->linename);
  bcopy (freqname, entry.line, 12);

  subscan = (int)(((headp->scan+0.005) - (int)headp->scan) * 100.0);

  strncpy(tmp, headp->backend, 8);
  tmp[8] = 0;                           /* String terminated backend name */

  if(verbose)
    st_report("Backend = %s, subscan %d", tmp, subscan );

  if(subscan < 20)
  {
    strcpy(telescope, "SMT-");
  }
  else
  if(subscan > 20 && subscan < 40)
  {
    strcpy(telescope, "HOT-");
  }
  else
  if(subscan > 40 && subscan < 60)
  {
    strcpy(telescope, "SKY-");
  }
  else
  if(subscan > 60 && subscan < 80)
  {
    strcpy(telescope, "CLD-");
  }
  else
    strcpy(telescope, "UNK-");

  sprintf(tbuf, "%s        ", cname+4);

  strcat(telescope, tbuf);

  bcopy((char *)telescope, (char *)&entry.telescope, 12);

  entry.date_obs = class_date(headp->utdate);
  entry.date_red = entry.date_obs;
  entry.offsets[0] = 0;                                      /* fill in below */
  entry.offsets[1] = 0;

  if(!strncmp(&headp->coordcd[4],"RADC", 4) || !strncmp(headp->coordcd,"APPRADC",7))
  {
    if(strncmp(headp->frame, "RA_DEC", 6) && strncmp(headp->frame, "      ", 6))
    {
      char cd[24], fr[24];

      strncpy( cd, headp->coordcd, 8 );
      cd[8] = 0;

      strncpy( fr, headp->frame, 8 );
      fr[8] = 0;

//      st_fail( "Inconsistent coordinate types: %s and %s", cd, fr );
//      st_fail( "Assume %4.4s is correct.", cd );
    }
    entry.offset_type = EQUAT;
  }
  else
  if(!strncmp(&headp->coordcd[0],"GALACTIC", 8))
  {
    entry.offset_type = GALAC;
  }
  else
  {
    st_fail("CLASS_WRITER: Unknown coordinate system %8.8s.", headp->coordcd);
    entry.offset_type = UNKNOWN;
  }

  entry.quality = 0;
  entry.scan = headp->scan;
  bzero((char *)entry.future, sizeof(entry.future));
                                                          /* fill scan header */
  scan_off = REC_LEN * (entry.addr-1);
  scan.entry = entry.obs;
  bcopy (SPECTRUM, scan.code, 4);
  scan.num_sec = 0;
  scan.sect_info[0] = SCAN_LEN + 1;
  scan.size = SCAN_LEN;
  scan.data_size = headp->datalen/4;
                                                   /* fill and write sections */
                                                           /* general section */
  scan.num_sec ++;
  add_class_section(scan.sect_info, scan.num_sec, GENERAL_SECT, GENERAL_LEN);
  scan.size += scan.sect_info[2*scan.num_sec-1];

  sect.general.ut_time = 3600*headp->ut;
  sect.general.lst_time = 3600*headp->lst;
  sect.general.azimuth = DEGREES2RADIANS(headp->az);
  sect.general.elevation = DEGREES2RADIANS(headp->el);
  sect.general.tau = headp->tauh2o + headp->tauo2;
  sect.general.tsys = headp->stsys;
  sect.general.integration = headp->inttime;

  sect_off = scan_off + (scan.sect_info[3*scan.num_sec-1]-1)*4;
  lseek(fd, sect_off, SEEK_SET);
  write_general(fd, &sect.general, scan.sect_info[2*scan.num_sec-1]*4);

                                                          /* position section */
  scan.num_sec ++;
  add_class_section(scan.sect_info, scan.num_sec, POSITION_SECT, POSITION_LEN);
  scan.size += scan.sect_info[2*scan.num_sec-1];

  bcopy(sourcename, sect.position.source, 12); 

  if(strcmp(headp->coordcd, "GALACTIC") == 0)
  {
    sect.position.epoch = 0;
  }
  else
  {
    if(strncmp(headp->coordcd, "APPRADC", 7) != 0 && 
                                atoi(headp->coordcd) != (int)headp->epoch)
    {
      char cd[24];

      strncpy( cd, headp->coordcd, 8 );
      cd[8] = 0;

      st_fail("Inconsistent coordinate epochs: %6.1f and %s", headp->epoch, cd);
      st_fail("Assume %6.1f is correct.", headp->epoch);
    } 
    sect.position.epoch = headp->epoch;
  }
                                                   /* Put in RA & Dec offsets */
  if(headp->openpar[8] != 9998.5 && headp->openpar[8] != 9999.5)
  {
    sect.position.lambda_off = 
     DEGREES2RADIANS((headp->xsource-headp->epocra)) * cos(sect.position.beta);

    sect.position.beta_off = 
     DEGREES2RADIANS((headp->ysource-headp->epocdec));
  }
  else                                 /* Az-El 5 pt.  Put in Az & El offsets */
  {
    sect.position.lambda_off = DEGREES2RADIANS(headp->az_offset);
    sect.position.beta_off   = DEGREES2RADIANS(headp->el_offset);
    entry.offset_type        = HORIZ;
  }

  if(entry.offset_type == GALAC)
  {
    sect.position.lambda = DEGREES2RADIANS(headp->gallong);       /* Galactic */
    sect.position.beta   = DEGREES2RADIANS(headp->gallat);

                                                /* Calculate Galactic offsets */
                    /* First get Galactic coordinates of the object + offsets */
    if(headp->openpar[8] != 9998.5 && headp->openpar[8] != 9999.5)
    {
      eqtogal(headp->xsource, headp->ysource, &glong_pos, &glat_pos);

      sect.position.lambda_off = 
          DEGREES2RADIANS(glong_pos - headp->gallong) * deg_cos(headp->gallat);

      sect.position.beta_off = 
          DEGREES2RADIANS(glat_pos - headp->gallat);
    }
  }
  else
  {
    sect.position.lambda = DEGREES2RADIANS(headp->epocra);        /* RA & Dec */
    sect.position.beta   = DEGREES2RADIANS(headp->epocdec);

                                                   /* Put in RA & Dec offsets */
    if(headp->openpar[8] != 9998.5 && headp->openpar[8] != 9999.5)
    {
      sect.position.lambda_off = 
      DEGREES2RADIANS((headp->xsource-headp->epocra)) * cos(sect.position.beta);

      sect.position.beta_off = 
       DEGREES2RADIANS((headp->ysource-headp->epocdec));
    }
    else                               /* Az-El 5 pt.  Put in Az & El offsets */
    {
      sect.position.lambda_off = DEGREES2RADIANS(headp->az_offset);
      sect.position.beta_off   = DEGREES2RADIANS(headp->el_offset);
      entry.offset_type        = HORIZ;
    }
  }

  entry.offsets[0] = sect.position.lambda_off;
  entry.offsets[1] = sect.position.beta_off;
  sect.position.projection = P_RADIO;

  sect_off = scan_off + (scan.sect_info[3*scan.num_sec-1]-1)*4;
  lseek (fd, sect_off, SEEK_SET);
  write_position(fd, &sect.position, scan.sect_info[2*scan.num_sec-1]*4);

                                                           /* spectro section */
  scan.num_sec ++;
  add_class_section (scan.sect_info, scan.num_sec, SPECTRO_SECT, SPECTRO_LEN);
  scan.size += scan.sect_info[2*scan.num_sec-1];

  bcopy(freqname, sect.spectro.line, 12);
  sect.spectro.rest = headp->restfreq;
  sect.spectro.nchan = scan.data_size;
  sect.spectro.reference = headp->refpt;
  sect.spectro.fres = headp->freqres;
  sect.spectro.foff = 0 ;                   /* what's proper interpretation ? */
  sect.spectro.vres = headp->deltax;

                                                 /* changed Apr 02, 2003: twf */
  if(headp->refpt_vel == 0.0)
  {
    sect.spectro.voff = headp->velocity;
  }
  else
  {
    sect.spectro.voff = headp->refpt_vel;
  }

  sect.spectro.bad = headp->badchv;
                                                 /*  sect.spectro.image = 0 ; */
  if(headp->sideband == 2.0)                                           /* USB */
  {
    sect.spectro.image = 
                 headp->restfreq - 2.0 * (headp->firstif + headp->openpar[0]);
  }

  if(headp->sideband == 1.0)                                           /* LSB */
  {
    sect.spectro.image = 
                 headp->restfreq + 2.0 * (headp->firstif - headp->openpar[0]);
  }
                                          /* ignore radio/optical definitions */
  if(!strncmp(&headp->veldef[4],"LSR", 3))
  {
    sect.spectro.velocity = VLSR;
  }
  else
  if(!strncmp(&headp->veldef[4], "HELO", 4)) 
  {
    sect.spectro.velocity = VHEL;
  }
  else
  if(!strncmp(&headp->veldef[4], "EART", 4)) 
  {
    sect.spectro.velocity = VEAR;
  }
  else
  if(!strncmp(&headp->veldef[4], "OBS", 3)) 
  {
    sect.spectro.velocity = VOBS;
  }
  else
  {
    st_fail("CLASS_WRITER: Unknown velocity reference %8.8s", headp->veldef);
    sect.spectro.velocity = VUNK;
  }

  sect.spectro.sky = headp->obsfreq;
  sect.spectro.vteles = headp->rvsys;                                    /* ? */
                                                        /* fix byte alignment */
  bcopy((char *)&sect.spectro.line,   sect_buf, 3*4);
  bcopy((char *)&sect.spectro.rest,  &sect_buf[3*4], 2*4);
  bcopy((char *)&sect.spectro.nchan, &sect_buf[(3+2)*4], 7*4);
  bcopy((char *)&sect.spectro.image, &sect_buf[(3+2+7)*4], 3*4); 
  bcopy((char *)&sect.spectro.sky,   &sect_buf[(3+2+7+3)*4], 3*4);

  sect_off = scan_off + (scan.sect_info[3*scan.num_sec-1]-1)*4;
  lseek (fd, sect_off, SEEK_SET);
  write(fd, (char *) &sect_buf, scan.sect_info[2*scan.num_sec-1]*4);

                                               /* frequency switching section */
  if(!strncmp(&headp->obsmode[4], "FS", 2))
  { 
   long nf;

    scan.num_sec ++;
    add_class_section(scan.sect_info, scan.num_sec, FSWITCH_SECT, FSWITCH_LEN);
    scan.size += scan.sect_info[2*scan.num_sec-1];

    sect.fswitch.nphase      = NPHASE;
    sect.fswitch.decalage[0] = headp->foffsig;
    sect.fswitch.decalage[1] = headp->foffref1;
    sect.fswitch.duree[0]    = headp->inttime/NPHASE;
    sect.fswitch.duree[1]    = headp->inttime/NPHASE;
    sect.fswitch.poids[0]    = 1;
    sect.fswitch.poids[1]    = -1;
    sect.fswitch.mode        = MOD_FREQ;
    sect.fswitch.lamoff[0]   = 0;
    sect.fswitch.betoff[0]   = 0;
    sect.fswitch.lamoff[1]   = 0;
    sect.fswitch.betoff[1]   = 0;
                             /* compress arrays & fix byte alignment */ 
    nf = NPHASE;
    bcopy((char *)&sect.fswitch.nphase, sect_buf, 4);
    bcopy((char *)&sect.fswitch.decalage, &sect_buf[4], nf*2*4 );
    bcopy((char *)&sect.fswitch.duree, &sect_buf[(1+nf*2)*4], nf*4);
    bcopy((char *)&sect.fswitch.poids, &sect_buf[(1+nf*3)*4], nf*4);
    bcopy((char *)&sect.fswitch.mode, &sect_buf[(1+nf*4)*4], 4);
    bcopy((char *)&sect.fswitch.lamoff, &sect_buf[(2+nf*4)*4], nf*4); 
    bcopy((char *)&sect.fswitch.betoff, &sect_buf[(2+nf*5)*4], nf*4); 

    sect_off = scan_off + (scan.sect_info[3*scan.num_sec-1]-1)*4;
    lseek(fd, sect_off, SEEK_SET);
    write(fd, (char *) &sect_buf, scan.sect_info[2*scan.num_sec-1]*4);
  }

                                                       /* calibration section */
  scan.num_sec ++;
  add_class_section(scan.sect_info, scan.num_sec, CALIBRATION_SECT, CALIBRATION_LEN);
  scan.size += scan.sect_info[2*scan.num_sec-1];

  sect.calib.beam_eff = headp->etal * headp->etafss;
  sect.calib.forw_eff = headp->etal;
  sect.calib.gain_im  = headp->gaini;
  sect.calib.h2omm    = headp->mmh2o;
  sect.calib.pamb     = headp->pressure;
  sect.calib.tamb     = headp->tamb + ZEROC;
  sect.calib.tatmsig  = headp->tatms;
  sect.calib.tchop    = headp->tchop;
  sect.calib.tcold    = headp->tcold;
  sect.calib.tausig   = headp->taus;     /* headp->tauh2o + headp->tauo2; */
  sect.calib.tauima   = headp->taui;
  sect.calib.tatmima  = 0; 
  sect.calib.trec     = headp->trx;
  sect.calib.mode     = 0;
  sect.calib.factor   = headp->tcal;

  if(!strncmp(headp->telescop,"NRAO 12M", 8))
    sect.calib.altitude = 1914;                                          /* m */
  else
    sect.calib.altitude = 3186; 

  sect.calib.counts[0] = headp->count[0];
  sect.calib.counts[1] = headp->count[1];
  sect.calib.counts[2] = headp->count[2];
  sect.calib.lamoff    = 0;
  sect.calib.betoff    = 0;

  sect_off = scan_off + (scan.sect_info[3*scan.num_sec-1]-1)*4;
  lseek(fd, sect_off, SEEK_SET);
  write(fd, (char *) &sect,scan.sect_info[2*scan.num_sec-1]*4);

                                                                /* write data */
  scan.data_addr =  scan.sect_info[3*scan.num_sec-1] + scan.sect_info[2*scan.num_sec-1];
  scan.size += scan.data_size;

  data_off = scan_off + (scan.data_addr-1)*4;
  lseek(fd, data_off, SEEK_SET);
  write(fd, (char *) datap, scan.data_size*4);

                         /* update the headers and indicies and close the file*/
  scan.size_rec = 1 + scan.size * 4 / REC_LEN;
  index.next_rec += scan.size_rec;

  lseek(fd, scan_off, SEEK_SET);
  write_scan(fd, &scan, sizeof (scan));

  ent_off = REC_LEN * (index.ext_addr [(scan.entry-1)/index.num_ent] +
    ((scan.entry-1) % index.num_ent)/ENT_REC - 1) + 
    ((scan.entry-1) % ENT_REC) * ENT_LEN;
  lseek(fd, ent_off, SEEK_SET);
  write_index(fd, &entry, sizeof (entry));

  lseek(fd, 0L, SEEK_SET);
  write_fileix(fd, &index, sizeof (index));
  close(fd);

  if(verbose)
  {
    st_report("Class: Obs type %8.8s, scan %7.2f written.", headp->obsmode, headp->scan);
  }

  return(1);
}


int class_init(file)  
char *file;
{
  int fd;
  struct file_index index;
  int len_index;

  if((fd = open(file, Oflag, Mode)) == -1)
  {
    st_fail("Can't open class file: %s", file);

    return(NOT_OK);
  }

  lseek(fd, 0L, SEEK_SET);
  len_index = read(fd, (char *) &index, sizeof (index));
  if(len_index == 0)
  {                                              /* file is empty; initialize */
    bcopy(EEEI_FILE, index.code, 4);
    index.next_rec = EXT_LEN / REC_LEN + 2 + 1;
    index.num_ent = ENT_EXT;
    index.num_ext = 1;
    index.next_ent = 1;
    index.ext_addr[0] = 3;
    bzero((char *) &index.ext_addr[1], sizeof (index.ext_addr)-4);
    write_fileix(fd, &index, sizeof (index));
    close(fd);

    return(NEW_OK);
  }
  else                                                         /* is file ok? */
  if(len_index == sizeof(index) && !bcmp(index.code, EEEI_FILE, 4))
  {
    close(fd);

    return(OLD_OK);
  }
  else                                                      /* file is not ok */
  {
    close(fd);
    st_fail("File %s isn't an EEEI CLASS file", file);

    return(NOT_OK);
  }
}



int write_index( fd, ix, len )
int fd;
struct index_ent *ix;
int len;
{
  struct index_ent xx;

  bcopy( ix, &xx, sizeof(struct index_ent));

  return( write( fd, &xx, len ));
}

int write_fileix( fd, ix, len )
int fd;
struct file_index *ix;
int len;
{
  struct file_index xx;

  bcopy( ix, &xx, sizeof(struct file_index));

  return( write( fd, &xx, len ));
}

int write_scan( fd, scan, len )
int fd;
struct scan_ent *scan;
int len;
{
  struct scan_ent ss;

  bcopy( scan, &ss, sizeof(struct scan_ent));

  return( write( fd, &ss, len ));
}

int write_general(fd, gen, len )
int fd;
struct general_descr *gen;
int len;
{
  struct general_descr gg;

  bcopy( gen, &gg, sizeof(struct general_descr));

  return( write( fd, &gg, len ));
}

int write_position(fd, gen, len )
int fd;
struct position_descr *gen;
int len;
{
  struct position_descr gg;

  bcopy( gen, &gg, sizeof(struct position_descr));

  return( write( fd, &gg, len ));
}