#include "templocator.h" #include "gondolageom.h" #include 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 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::getAzimutSST(int sampleNum) { findGeomFromGC(sampleNum); double az = azimBolGC; if (az>360) az -= 360; if (az<0) az += 360; return az; } double TempLocator::getElvSST(int /*sampleNum*/) { return GondolaGeom::elevSST0; } double TempLocator::getAlphaSST(int sampleNum) { double elv = getElvSST(sampleNum); double az = getAzimutSST(sampleNum); 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::getDeltaSST(int sampleNum) { double elv = getElvSST(sampleNum); double az = getAzimutSST(sampleNum); double ha; //double ra; double 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 dec; } 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; double 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 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