source: Sophya/trunk/Poubelle/archTOI.old/templocator.cc@ 534

#include "templocator.h"
#include <math.h>
extern "C" {
#include "aa_hadec.h"
}

#include "fitsio.h"
#include "plgalcross.h"
#include "archparam.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

TempLocator tempLocator;

TempLocator::TempLocator()
{
  lon = lat = ts = 0;
  raZ  = decZ = -99999;
  xSampleNum = -99999;
  fitsfile* fptr;
  int status=0;
  fits_open_file(&fptr, "samplenum_gal_cross.fits", READONLY, &status);
  int simple, bitpix, naxis;
  long naxes;
  long pcount, gcount;
  int extend;
  fits_read_imghdr(fptr, 1, &simple, &bitpix, &naxis, &naxes, &pcount, &gcount, &extend, &status);
  nGalCross = naxes;
  crossings = new long[nGalCross];
  int anynul;
  fits_read_img_lng(fptr, 0, 1, nGalCross, 0, crossings, &anynul, &status);
  fits_close_file(fptr, &status);
  fits_report_error(stderr, status);  /* print out any error messages */
}
           
           
void TempLocator::setEarthPos(double lon, double lat) {
  if (this->lon == lon && this->lat == lat) return;
  this->lon = lon;
  this->lat = lat;
  raZ  = decZ = -99999; xSampleNum = -99999;
}

void TempLocator::setTSid(double ts) {
  if (this->ts == ts) return;
  this->ts  = ts;
  raZ  = decZ = -99999; xSampleNum = -99999;
}
  
void   TempLocator::ComputeZenith() {
  double ha;
  aa_hadec (lat * M_PI/180, .5 * M_PI, 0, &ha, &decZ);
  raZ = - (ha * 180. / M_PI / 15) + (ts/3600.);
  if (raZ>24) raZ -= 24;
  if (raZ<0)  raZ += 24;
  decZ = decZ * 180. / M_PI;
}
  
double TempLocator::getAlphaZenith() {
  if (raZ < -100) ComputeZenith();
  return raZ;
}

double TempLocator::getDeltaZenith() {
  if (decZ < -100) ComputeZenith();
  return decZ;
}

#define altbolo1 41.5

void TempLocator::findGeomFromGC(int sampleNum) // pour le bolo qui voit les xing
{
  if (sampleNum == xSampleNum) return;
  if (decZ < -100) ComputeZenith();
  
  azimBolGC = -9999;
  
  // On trouve les croisements juste avant et juste apres notre sampleNum
  for (icross=0; icross<nGalCross; icross++) {
    if (crossings[icross] > sampleNum) break;
  }
  if (icross == 0 || icross >= nGalCross) return;
  
  // On trouve l'azimut du croisement principal pour notre position actuelle
  double alpG = 12. + 51./60. + 30./3600.;
  double delG = 27. + 07./60. + 42./3600.;
  double azCr1, azCr2;
  int rc = PlGalCross(ts/3600., lat, (90. - altbolo1), alpG, delG, azCr1, azCr2);
  if (rc != 0) return; // pas deux points d'intersection
  
  // Il faut determiner le croisement principal, ie le plus proche
  // du centre galactique. Pendant le vol de Trapani, c'etait celui
  // le plus proche de 220° d'azimut.
  
  double azCross = azCr1;
  if (fabs(azCr2-220) < fabs(azCr1-220)) azCross = azCr2;
  
  rotSpeed = 360./(crossings[icross] - crossings[icross-1]); // °/sample
  
  azimBolGC = azCross - (sampleNum - (crossings[icross-1]+12))*rotSpeed;
  // azimut bolo 1 from crossing, for sampleNum
  if (azimBolGC > 360) azimBolGC -= 360;
}

double TempLocator::getRotSpeed(int sampleNum) {
  findGeomFromGC(sampleNum);
  return rotSpeed / archParam.acq.perEch;
}

int  TempLocator::getCrossSamples(int sampleNum, int& SN1, int& SN2) {
  findGeomFromGC(sampleNum);
  if (icross == 0 || icross >= nGalCross) return -1;
  SN1 = crossings[icross-1];
  SN2 = crossings[icross];
  return 0;
}


void TempLocator::getAltAzBolo(int sampleNum, int ibolo, double& elv, double& az) {
  findGeomFromGC(sampleNum);
  double delElv = 0;
  double delAz = 0;  // relative to ch1-bolo1
  elv = -99999;
  az  = -99999;
  switch (ibolo) {
    case 11:
      delElv =   0;
      delAz  =   0;
      break;
    case 8:
      delElv =   0.78 * sqrt(3.)/2.;
      delAz  =   0.78 * 1./2.;
      break;
    case 13:
      delElv = - 0.78 * sqrt(3.)/2.;
      delAz  =   0.78 * 1./2.;
      break;
    case 9:
      delElv =   0.78 * sqrt(3.)/2.;
      delAz  = - 0.78 * 1./2.;
      break;
    case 4:
      delElv =   0.;
      delAz  =   0.78;
      break;
    case 15:
      delElv =   0.;
      delAz  = - 0.78;
      break;
    default:
      return;
  }
  delAz /= cos(41 * M_PI/180);
  elv = altbolo1 + delElv;
  az  = azimBolGC + delAz;
  if (az>360) az -= 360;
  if (az<0)   az += 360;
  return;
}

double TempLocator::getAzimutBolo(int sampleNum, int ibolo) {
  double elv, az;
  getAltAzBolo(sampleNum, ibolo, elv, az);
  return az;
}

double TempLocator::getElvBolo(int sampleNum, int ibolo) {
  double elv, az;
  getAltAzBolo(sampleNum, ibolo, elv, az);
  return elv;
}

double TempLocator::getAlphaBolo(int sampleNum, int ibolo) {
  double elv, az;
  getAltAzBolo(sampleNum, ibolo, elv, az);
  double ha;
  double ra,dec;
  aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
  ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
  while (ra>24) ra -= 24;
  while (ra<0) ra += 24;
  dec = dec * 180. / M_PI;
  return ra;
}

double TempLocator::getDeltaBolo(int sampleNum, int ibolo) {
  double elv, az;
  getAltAzBolo(sampleNum, ibolo, elv, az);
  double ha;
  double ra,dec;
  aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
  ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
  dec = dec * 180. / M_PI;
  return dec;
}

double TempLocator::getAlphaCenter(int sampleNum) {
  return getAlphaBolo(sampleNum, 11);
}

double TempLocator::getDeltaCenter(int sampleNum) {
  return getDeltaBolo(sampleNum, 11);
}

double TempLocator::getAzimutCenter(int sampleNum) {
  return getAzimutBolo(sampleNum, 11);
}

double TempLocator::getElvCenter(int sampleNum) {
  return getElvBolo(sampleNum, 11);
}




                        //       ;
                        //       ;                            /\         4     2
                        //       ;                  elevation ||                       positive scanning, clockwise
                        //       ;                            ||       6    1    5      ---------->
                        //       ;                            ||                      == positive azimut
                        //       ;                            ||         x     3
            //                                         |----|      = 0.78 deg / cos(elev)
            // bol 1 = 41° elevation