/*+ e c h T e l M d l . c * Module name: echTelMdl.c * Function: Echidna-FMOS Telescope Model library. * Description: This library provides the FMOS/Echida Telescope Model Implementation. The library implements features to provide for transformation of sky positions in RA/DEC coordinates to Echidna positioner coordinates. This file is used by both the instrument software and the Echidna configure (fibre allocation) software. The following classes are declared in the header echTelMdl, namespace echidna telMdlException => The exception class thrown by these classes. A sub-class of EchidnaException. telMdlPM_PARS => A simple struct used for passing around proper motion details. The contructor tasks the ra, dec, parallax and radialVel proper motion values in that order, each value defaulting to zero. telMdlFlexCorPars => Parameters involved in the correction of telescope tube flexure. Details TBD. telMdlPARS => A structure containing an array of doubles of size telMdlNUM_PARS. This is used to passing the telescope model parameters around for for allowing access to the indivual values (required by the calibration software). Access to the indivual parameters is supported using the .p item. telMdl => A class which provides access to the telescope model. Details below. telMdlMean => A sub class of telMdl which allows for conversion of mean positions rather the just apparent positions. The following is an introduction to the telMdl public methods. Two constructors a provided. On takes a telMdlPARS item as a argument and initialises the model from those parameters. The other loads the model from a file - with the filename and directory specfications defaulting appropiately. Before a conversion can be done, the model must be initialised for the current astmospheric condition and the field centre. This is done using the CvtInit() method. The method App2Pos() will convert an apparent position on the sky to a echidna x/y positions. The method App2Tan() will convert an apparent position on the sky to tangent plane coorindates. The Pos2App() method will convert echidna x/y positions to an approximate sky position. The method QuickAppToObs() does a quick conversion from apparent to observed coorindates. The RA and Dec passed to CvtInit() can be retreived using GetCenRA() and getCenDec() (particually usefull when the telMdlMean sub-class is used, as this allows access to the apparent field center position). The method Get() allows the telescope model parameters to be fetched, whilst Set() allows them to be set. The method Normalize() is used to remove the extra scale, rotation and non-perpendicularity in a model by normalizing them into the base linear and distortion model. The Micros2Arssec() method is used to conversion microns on the field plate to an approimaton of arc-seconds on the sky. The method Write() will write the model to a file. The directory name and be optinally specified. * Language: C++ * Namespace:echidna * Support: Tony Farrell, AAO *- * $Id: ACMM:echTelMdl/src/echTelMdl.cc,v 2.4 29-Jan-2008 01:07:08+11 tjf $ * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz C version. 12-Jun-2003 - SS - set WAVE_MAX to 1.8 micron so I can test with some realistic IR values. 10-Jun-2004 - SS - Make Pos2App work (remove hardcoded constants, etc.); Fix corrector distortion model; Created Mean2Pos(). 17-Jun-2004 - SS - Added Pos2AppIter() & Pos2MeanIter(). 18-Jun-2004 - SS - Add optional parameter $equinox to Pos2MeanIter() & Pos2Mean(); Make member-functions const where possible. 24-Jan-2006 - SS - Fix XY alignment of flexure compensation. 05-Sep-2006 - SS - bugfix: usage of telMdlPAR_Y_SCALE & telMdlPAR_Y_SHEAR was incorrect - they were swapped everywhere so the code worked but shear actually referred to scale & vice-versa. 10-Oct-2006 - SS - Use $TEL_MDL_FOV_WARN for threshold value. 08-Feb-2007 - SS - Created $USE_OLD_DIST_MDL for debugging telescope model with OzPos data. * Copyright (c) Anglo-Australian Telescope Board, 2007. Not to be used for commercial purposes without AATB permission. */ static const char *rcsId="@(#) $Id: ACMM:echTelMdl/src/echTelMdl.cc,v 2.4 29-Jan-2008 01:07:08+11 tjf $"; static void *use_rcsId = ((void)&use_rcsId,(void *) &rcsId); // Include files. #include #include #include #include #include extern "C" { #include "slalib.h" #include "slamac.h" } #include "echTelMdl.h" #include "log.hh" #include "envVars.hh" // getStringEV() getFloatEV() #include "macros.hh" #include #include #include using std::cout; using std::cerr; using std::endl; /* * Some macros to ease access to exceptions. THROW just throws with * a given mesage. EXCEPTION_STREAM creates an exception stream object * which is thrown when the destructor is invoked. */ #define THROW(_message...) throw telMdlException(__FILE__,__LINE__,##_message) #define EXCEPTION_STREAM(_obj) echidna::exceptionStream _obj(__FILE__,__LINE__) /* * Parameter ranges for the parameters to telMdlCvtInit(). */ const double TEMP_MIN = 200.0; /* Min temperature we expect to see -73 C */ const double TEMP_MAX = 330.0; /* Max temperature we expect to see +57 C */ const double PRES_MIN = 200.0; /* Min pressure we expect to see - millibars */ const double PRES_MAX = 1000.0;/* Max pressure we expect to see - millibars */ const double HUMI_MIN = 0.0; /* Min humidity we expect to see - relative (0%) */ const double HUMI_MAX = 1.0; /* Max humidity we expect to see - relative (100%) */ const double WAVE_MIN = 0.3; // Min wavelength we expect to see - µm (300 nm) const double WAVE_MAX = 1.8; // Max wavelength we expect to see - µm (1800 nm) /* * Maximum X/Y value - 1000000 microns, 1 meter, well outside the plate. */ const unsigned X_Y_MAX = 1000000; namespace echidna { const char * const telMdl::MODEL_FILE = "telescopeModel.txt"; const bool bDebugTM = #ifdef TESTCODE getBoolEV("DEBUG_TELESCOPE_MODEL", false); #else false; #endif const double TRP_LSP_RATE = 0.0065;// Tropospheric lapse rate - see slaAop etc. /* * Class Private Function name: telMdl::Dump * Description: Debugging function - dumps the model to stderr. * History: 20-Feb-2003 - TJF - Original version */ void telMdl::Dump() const { fprintf(stderr,"Dump of telescope model info, model address = %p\n", this); fprintf(stderr,"Distortion = "); for (unsigned i = 0 ; i < telMdlNUM_DIST_PARS ; ++i) fprintf(stderr,"%g ", distPars[i]); fprintf(stderr,"\n"); fprintf(stderr,"Linear Fwd = %g %g %g %g %g %g\n", linPars.coeffs[0], linPars.coeffs[1], linPars.coeffs[2], linPars.coeffs[3], linPars.coeffs[4], linPars.coeffs[5]); fprintf(stderr,"Linear Inv = %g %g %g %g %g %g\n", linPars.invCoeffs[0], linPars.invCoeffs[1], linPars.invCoeffs[2], linPars.invCoeffs[3], linPars.invCoeffs[4], linPars.invCoeffs[5]); fprintf(stderr,"Scale = %g, Rotation = %g (degrees), Non-perp = %g(degrees)\n", linPars.scale, linPars.rotation, linPars.nonperp); fprintf(stderr,"XY cenAopars = 0(%g) 1(%g) 2(%g) 3(%g) 4(%g)\n", cenAoprms[0], cenAoprms[1], cenAoprms[2], cenAoprms[3], cenAoprms[4]); fprintf(stderr," 5(%g) 6(%g) 7(%g) 8(%g) 9(%g)\n", cenAoprms[5], cenAoprms[6], cenAoprms[7], cenAoprms[8], cenAoprms[9]); fprintf(stderr," 10(%g) 11(%g) 12(%g) 13(%g)\n", cenAoprms[10], cenAoprms[11], cenAoprms[12], cenAoprms[13]); fprintf(stderr,"XY obsAopars = 0(%g) 1(%g) 2(%g) 3(%g) 4(%g)\n", obsAoprms[0], obsAoprms[1], obsAoprms[2], obsAoprms[3], obsAoprms[4]); fprintf(stderr," 5(%g) 6(%g) 7(%g) 8(%g) 9(%g)\n", obsAoprms[5], obsAoprms[6], obsAoprms[7], obsAoprms[8], obsAoprms[9]); fprintf(stderr," 10(%g) 11(%g) 12(%g) 13(%g)\n", obsAoprms[10], obsAoprms[11], obsAoprms[12], obsAoprms[13]); fprintf(stderr,"XY Cen lambda = %g\n", cenLambda); fprintf(stderr,"XY Obs lambda = %g\n", obsLambda); fprintf(stderr,"XY ypa = %g\n", ypa); fprintf(stderr,"---------------------------\n"); } /* * Class Private Function name: telMdl::OpenModelFile * Description: Open a file containing the Echidna transformation details. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) filename (const std::string *) Base name of file to open. If NULL then use the standard name. (>) dirOverride (const std::string *) If non-null, a directory to search for the distortion details file. (>) dirFallback (const std::string *) If non-null, a fall back directory to search for the distortion details file. (!) file (ifstream &) The stream in which the file will be opened. (<) actfilename (std::string *) The actual name of the file opened will be returned here. * History: 20-Feb-2003 - TJF - Original version */ void telMdl::OpenModelFile ( const std::string * filename, const std::string * const dirOverride, const std::string * const dirFallback, std::ifstream & file, std::string * const actfilename) const { std::string default_name(MODEL_FILE); if (!filename) filename = &default_name; /* * Try to open the file in dirOverride if one was specified */ if (dirOverride != 0) { *actfilename = *dirOverride + "/" + *filename; file.open(actfilename->c_str()); if (file) return; } /* If not yet found try to open in default directory */ *actfilename = *filename; file.open(actfilename->c_str()); if (file) return; /* If still not found try to open in dirFallback */ if (dirFallback) { *actfilename = *dirFallback + "/" + *filename; file.open(actfilename->c_str()); if (file) return; } THROW("Could not find telescope model file '%s'", filename->c_str()); } /* *+ t e l M d l : : t e l M d l * Function name: telMdl::telMdl * Function: Construct a telMdl object using supplied parameters. * Description: This method initializes access to the Echidna transformation details from the supplied parameters. * Language: C++ * Namespace:echidna * Call: telMdl::teMdl(pars) * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) pars (const telMdlPARS &) Contains the parameters of the telescope model to use. * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ telMdl::telMdl(const telMdlPARS &pars) { CreateModel(pars); } /* *+ t e l M d l : : t e l M d l * Function name: telMdl::telMdl * Function: Construct a telMdl object, inputting parameters from a file. * Description: This method initializes access to the Echidna transformation details from calibration model files (named EchidnaModel.txt) * Language:C++ * Namespace:echidna * Call: telMdl::teMdl(filename=0, dirOverride=0, dirFallback=0) * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) filename (const std::string *) If non-null, the base name of file to open. If null, then the standard name is used. (>) dirOverride (const char *) If non-null, a directory to search for the distortion details file. (>) dirFallback (const char *) If non-null, a fall back directory to search for the distortion details file. * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. 10-Dec-2003 - SS - reading of config file wasn't working with g++-3.2 */ telMdl::telMdl ( const std::string * const filename, const std::string * const dirOverride, const std::string * const dirFallback) { /* * Something is wrong if we have less then 6 distortion parameters */ assert(telMdlNUM_DIST_PARS >= 6); /* * Need a string for the file name. */ std::string filebuf; /* * The input stream we will use. */ std::ifstream file; /* * Read the telescope model information from the file. * * First open the file. */ OpenModelFile( filename, dirOverride, dirFallback, file, &filebuf); std::cout << "Reading telescope model parameters from " << filebuf << std::endl; telMdlPARS pars; /* * Read the contents of the model (need to improve error checking) */ std::string s; for (unsigned i = 0; i < telMdlNUM_PARS ; ++i) { if (file.eof()) THROW("Error reading parameter %d from telescope model file %s", i+1 , filebuf.c_str()); getline(file, s); pars.p[i] = atof(s.c_str()); } /* * Create the model. */ CreateModel(pars); } /* * Class Private Function name: telMdl::CreateModel * Description: This method initialises access to the Echidna transformation details. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) pars (const telMdlPARS &) Contains the parameters of the telescope model to use. * History: 20-Feb-2003 - TJF - Original version 05-Sep-2006 - SS - bugfix: usage of telMdlPAR_Y_SCALE & telMdlPAR_Y_SHEAR was incorrect 25-Jan-2008 - TJF - Use Scott's FEQUAL macros to check floating point equality. */ void telMdl::CreateModel ( const telMdlPARS ¶meters) { double a,b,c,d,e,f; double det; const double *pars = ¶meters.p[0]; // Easy access to parameters. /* Set up linear coefficients */ a = pars[telMdlPAR_X_ORIGIN]; b = pars[telMdlPAR_X_SCALE]; c = pars[telMdlPAR_X_SHEAR]; d = pars[telMdlPAR_Y_ORIGIN]; e = pars[telMdlPAR_Y_SHEAR]; f = pars[telMdlPAR_Y_SCALE]; linPars.coeffs[0] = a; linPars.coeffs[1] = b; linPars.coeffs[2] = c; linPars.coeffs[3] = d; linPars.coeffs[4] = e; linPars.coeffs[5] = f; /* Calculate coefficients for inverse transformation */ det = b * f - c * e; if (FEQUAL9(det,0)) // Check to 9 decimal places THROW("Invalid linear transformation - no inverse"); linPars.invCoeffs[0] = ( c * d - a * f ) / det; linPars.invCoeffs[1] = f/det; linPars.invCoeffs[2] = -1 * c/det; linPars.invCoeffs[3] = ( a * e - b * d ) / det; linPars.invCoeffs[4] = -1 * e / det; linPars.invCoeffs[5] = b / det; /* Rotation, scale and non-perpendicularity for the linear model */ linPars.rotation = pars[telMdlPAR_ROT]; linPars.scale = pars[telMdlPAR_SCALE]; linPars.nonperp = pars[telMdlPAR_NONPERP]; /* The rest - considered to be distortion */ unsigned i,j; for ( i = 0, j = telMdlPAR_NW_DIST_FIRST ; i < telMdlNUM_DIST_PARS ; ++i, ++j) distPars[i] = pars[j]; } /* *+ t e l M d l : : G e t * Function name: telMdl::Get * Function: Return the parameters of a model * Description: This method returns an array containing the parameters of a model in the same format as can be passed to a telMdl::telMld() constructor. * Language: C++ * Namespace:echidna * Call: (telMdlPARS) = telMdl::Get() const * Returned value: The telescope model parameters. * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. 19-Dec-2005 - SS - Unhardcode constants - use telMdlPAR_* names. 05-Sep-2006 - SS - bugfix: usage of telMdlPAR_Y_SCALE & telMdlPAR_Y_SHEAR was incorrect */ telMdlPARS telMdl::Get() const { unsigned i,j; telMdlPARS p; double *pars = &p.p[0]; /* Zero all values */ for (i = 0 ; i < telMdlNUM_PARS ; ++i) pars[i] = 0; /* * Get the linear model values. */ pars[telMdlPAR_X_ORIGIN] = linPars.coeffs[0]; pars[telMdlPAR_X_SCALE] = linPars.coeffs[1]; pars[telMdlPAR_X_SHEAR] = linPars.coeffs[2]; pars[telMdlPAR_Y_ORIGIN] = linPars.coeffs[3]; pars[telMdlPAR_Y_SHEAR] = linPars.coeffs[4]; pars[telMdlPAR_Y_SCALE] = linPars.coeffs[5]; pars[telMdlPAR_ROT] = linPars.rotation; pars[telMdlPAR_SCALE] = linPars.scale; pars[telMdlPAR_NONPERP] = linPars.nonperp; /* * Fetch the distortion model values. */ for (i = 0, j = telMdlPAR_NW_DIST_FIRST ; i < telMdlNUM_DIST_PARS ; ++i, ++j) pars[j] = distPars[i]; return p; } /* *+ t e l M d l : : S e t * Function name: telMdl::Set * Function: Set the parameters of a model * Description: This method is given an array containing the parameters of a model and changes the parameters of the model. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::Set(pars) * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) pars (const telMdlPARS &) Array containing model parameters. * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ extern void telMdl::Set( const telMdlPARS &p) { // Just call CreateModel CreateModel(p); } /* *+ t e l M d l : : W r i t e * Function name: telMdl::Write * Function: Write model to a file. * Description: This method writes the model to file of suitable format for picking up by the telMdl constructor which reads from files. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::Write(dirname = 0) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) dirname (const string *) Name of directory into which file will be written. If null, the current directory is used. * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ void telMdl::Write( const std::string * const dirname) const { std::string filename; /* * Get the full file name. */ if (dirname != 0) filename = *dirname + "/" + MODEL_FILE; else filename = MODEL_FILE; /* * Create file. */ std::ofstream file(filename.c_str()); if (!file) THROW("Failed to create file %s", filename.c_str()); /* * Write parameters. */ telMdlPARS pars = Get(); for (unsigned i = 0; i < telMdlNUM_PARS; ++i) { file << pars.p[i] << '\n'; } } /* * Class Private Function name: telMdl::ValidatePars * Description: Validate the values of parameters passed to telMdl::CvtInit(). * History: 20-Feb-2003 - TJF - Original version 10-Dec-2003 - SS - split function into 2 - new checkObservable(). */ void telMdl::ValidatePars ( double mjd, double temp, double press, double humid, double cenWave, double obsWave, double rotation, double ra, double dec) const { if ((ra < 0)||(ra > D2PI)) { THROW( "Invalid field center RA value (%.6f radians)",ra); } else if ((dec < -DPI)||(dec > DPI)) { THROW("Invalid field center Dec value (%.6f radians)", dec); } else if ((rotation < -2*M_PI)||(rotation > 2*M_PI)) { THROW("The rotator position angle of %g radians is outside the expected range", rotation); } /* * Validate passed parameters. */ else if ((temp < TEMP_MIN)||(temp > TEMP_MAX)) { THROW("The tempreature value of %g kelvin is outside the expected range of %g to %g", temp, TEMP_MIN, TEMP_MAX); } else if ((press < PRES_MIN)||(press > PRES_MAX)) { THROW("The pressure value of %g miliibars is outside the expected range of %g to %g", press, PRES_MIN, PRES_MAX); } else if ((humid < HUMI_MIN)||(humid > HUMI_MAX)) { THROW("The relative humidity value of %g is outside the expected range of %g to %g", humid, HUMI_MIN, HUMI_MAX); } else if ((cenWave < WAVE_MIN)||(cenWave > WAVE_MAX)) { THROW("The centre wavelength value of %.3g µm is outside the expected range of %.3g to %.3g", cenWave, WAVE_MIN, WAVE_MAX); } else if ((obsWave < WAVE_MIN)||(obsWave > WAVE_MAX)) { THROW("The observation wavelength value of %.3g µm is outside the expected range of %.3g to %.3g", obsWave, WAVE_MIN, WAVE_MAX); } (void) mjd; // unused parameter - pacify compiler. } /* * Internal Function name: telMdl::checkObservable * Description: Check that the field is observable. ie. above the horizon. * History: 10-Dec-2003 - SS - Initial version. 24-Jan-2008 - TJF - Remove old comment about ozd > 1.57 check not working. It is working now. */ void telMdl::checkObservable ( double mjd, double ra, double dec) { /* * If the field is essentially * unobservable because of an unrealistic mjd value, then slaAopqk() * will take forever doing the refraction calculations. So we do an * approximate calculation first, disabling refraction by setting the * pressure parameter to zero. If the calculated zenith distance is * less than 90 degrees (~ 1.57 rad) we know the full calculation is * feasible and repeat it properly. * * Create a copy of the paramters and set the pressure to zero. */ double modCenAoprms[numAoPars]; memcpy(modCenAoprms, cenAoprms, sizeof(modCenAoprms)); slaAoppat(mjd, modCenAoprms); modCenAoprms[6] = 0.0; /* * Do calculation */ double oaz,ozd,oha, ocdec, ocra; slaAopqk(ra,dec,modCenAoprms,&oaz,&ozd,&oha,&ocdec,&ocra); if (bDebugTM) cout << "Zenith Distance=" << (ozd*180.0/M_PI) << " degrees." << endl; /* * Catch and report error. */ if (ozd > 1.57) { EXCEPTION_STREAM(e); // Exception thrown on exit from block. char ra_sign; char dec_sign; char time_sign; int hmsf[4]; int dmsf[4]; int ymdf[4]; int thmsf[4]; int j; e << "Illegal date of observation, field details follow" << std::endl; e << "The MJD of " << std::fixed << std::showpoint << std::setprecision(6) << mjd << " leads to invalid zenith distance of " << std::setprecision(3) << ozd << " rad (" << std::setprecision(0) << (ozd*180.0/M_PI) << " deg)" << std::endl; slaDr2tf(2, ra, &ra_sign, hmsf); slaDr2af(2, dec, &dec_sign, dmsf); slaDjcal(2, mjd, ymdf, &j); slaCd2tf(0, mjd - (int)mjd, &time_sign, thmsf); e << "Field Centre (app) RA = " << std::setprecision(6) << ra << "(" << ra_sign << std::setw(2) << hmsf[0] << ":" << hmsf[1] << ":" << hmsf[2] << "." << hmsf[3] << ") Dec = " << dec << "(" << dec_sign << dmsf[0] << ":" << dmsf[1] << ":" << dmsf[2] << "." << dmsf[3] << ")" << std::endl; e << "Observation Time:MJD = " << std::setprecision(8) << mjd << ", UT = " << std::setw(4) << ymdf[0] << "/" << std::setw(2) << ymdf[1] << "/" << ymdf[2] << " " << thmsf[0] << ":" << thmsf[1] << ":" << thmsf[2]; // Exception now thrown in destructor } } /* * Class Private Function name: telMdl::LogCvtInit * Description: Log the values of the parameters passed to CvtInit. * History: 20-Feb-2003 - TJF - Original version */ void telMdl::LogCvtInit ( double mjd, double temp, double press, double humid, double cenWave, double obsWave, double rotation, double ra, double dec) const { char buffer[200]; char time_sign; int ymdf[4]; int thmsf[4]; int j; slaDjcal(2, mjd, ymdf, &j); slaCd2tf(0, mjd - (int)mjd, &time_sign, thmsf); sprintf(buffer, "Telescope model initialised at mjd=%.6g (%.4d/%.2d/%.2d %.2d:%.2d:%.2d)", mjd, ymdf[0], ymdf[1], ymdf[2], thmsf[0], thmsf[1], thmsf[2]); llog << buffer << std::endl; /* */ char ra_sign; char dec_sign; int hmsf[4]; int dmsf[4]; slaDr2tf(2, ra, &ra_sign, hmsf); slaDr2af(2, dec, &dec_sign, dmsf); sprintf(buffer, "Supplied Field Centre (app) RA = %.6f(%c%.2d:%.2d:%.2d.%d) Dec = %.6f(%c%.2d:%.2d:%.2d.%d)", ra, ra_sign, hmsf[0], hmsf[1], hmsf[2], hmsf[3], dec, dec_sign, dmsf[0], dmsf[1], dmsf[2], dmsf[3]); llog << buffer << std::endl; sprintf(buffer, "Telescope model initialised with cwave=%.3f nm, owave=%.3f nm temp=%g K, pres=%g mb, humid=%g", cenWave, obsWave, temp, press, humid); llog << buffer << std::endl; sprintf(buffer,"Position angle = %.6f", rotation); llog << buffer << std::endl; } /* *+ t e l M d l : : C v t I n i t * Function name: telMdl::CvtInit * Function: Initialize for conversions at a specified time, place etc. * Description: This method initializes the class for conversions at a particular time, wavelength and telescope model details and apparent place. This allows a sequence of transformations for one field to be done efficiently. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::CvtInit(mjd, dut, temp, pres, humid, cenWave, obsWave, rotation, flexCorPars, ra, dec) * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) mjd (double) UTC date and time as a modified Julian date. (>) dut (double) Delta UT (UT1 - UTC) seconds. (>) temp (double) Atmospheric temperature (K). (>) press (double) Atmospheric pressure (mB). (>) humid (double) Atmospheric humidity (0 to 1.0). (>) cenWave (double) Wavelength telescope is pointed with (microns) (>) obsWave (double) Wavelength of observation (micons). (>) rotatorPA (double) Rotator position angle (radians). (>) flexCorPars (const telMdlFlexCorPars &) flex correction parameters. (>) ra (double) Apparent RA of observation (radians). (>) dec (double) Apparent Dec of observation (radians) * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. 10-Dec-2003 - SS - use $OBSERVATORY_NAME environment variable. bugfix: checkObservable() was being called too early. 24-Jan-2008 - TJF - Longitude, latitude and height of observatory are now stored in class variables. Calculate rotator setting. */ void telMdl::CvtInit ( double mjd, double temp, double press, double humid, double cenWave, double obsWave, double rotatorPA, const telMdlFlexCorPars &flexCorPars, double ra, double dec) { if (bDebugTM) LogCvtInit(mjd, temp, press, humid, cenWave, obsWave, rotatorPA, ra, dec); char name[80]; /* * Validate Parameters. */ ValidatePars(mjd, temp, press, humid, cenWave, obsWave, rotatorPA, ra, dec); if (humid < 0.001) { /* * Just in case slaLib has trouble with humidty of zero, which is not * impossible after rounding etc, set a humidity of zero * to 0.001. */ humid = 0.001; } /* * Save these field center values. */ dateMjd = mjd; cenRa = ra; cenDec = dec; dut = 0; /* Normally don't bother setting this */ /* * Get observatory parameters */ static const char *obsName = getStringEV("OBSERVATORY_NAME", "SUBARU"); slaObs(-1,const_cast(obsName),name,&longit_,&lat_,&height_); if (strcmp(name,"?") == 0) THROW("\"%s\" telescope not recognised by the version of SLA being used", obsName); /* * Determine apparent to observed place parameters */ slaAoppa(mjd,dut,-longit_,lat_,height_,0.0,0.0,temp,press,humid, cenWave, TRP_LSP_RATE, cenAoprms); slaAoppa(mjd,dut,-longit_,lat_,height_,0.0,0.0,temp,press,humid, obsWave, TRP_LSP_RATE, obsAoprms); /* * Save the wavelengths. */ cenLambda = cenWave; obsLambda = obsWave; /* * params[0] is the argus orientation. We add the default * rotatorPA to this to get the rotation to be applied. */ ypa = -rotatorPA; /* * Determine the zenith distance - required in the distortion calculation. */ double aob, hob, dob, rob; QuickAppToObs(ra, dec, &aob, &zenithDistance, &hob, &dob, &rob); /* * Detemrine the rotator setting at this point (used by GetDistCenter(). * We just add the parallactic angle to the position angle */ rotatorSetting_ = slaPa(hob, dob, lat_) + rotatorPA; /* * Save this. (not yet used (24-Jan-2008). */ flexCorrectionPars = flexCorPars; /* * Check we can observe this field. */ checkObservable(mjd, ra, dec); } /* * Class Private Function name: telMdl::ValidateApp2Tan * Description: Validate an apparent RA/Dec sky position. * History: 20-Feb-2003 - TJF - Original version */ void telMdl::ValidateApp2Tan( double ra, double dec) const { char mesbuf[200]; int bad=0; EXCEPTION_STREAM(e); /* Exception thrown on exit from block. unless dontThrow invoked */ /* * Default settings for exception stream and first part of message. */ e << std::fixed << std::showpoint << std::setprecision(6); e << "Invalid arguments range in call to telMdl::ValidateApp2Tan"; if ((ra > D2PI) || (ra < 0.0)) { bad++; e << "Invalid RA value (" << ra << " radians)" << std::endl; } if ((dec > DPI) || (dec < (-DPI))) { bad++; e << "Invalid Dec value (" << dec << " radians)" << std::endl; } if (bad) { char ra_sign; char dec_sign; int hmsf[4]; int dmsf[4]; slaDr2tf(2, ra, &ra_sign, hmsf); slaDr2af(2, dec, &dec_sign, dmsf); //#warning "We should have a user type with RA/Dec for which the formating is done" // SS: There is an 'RaDec' class in the echGenericModules module. // Currently the 'operator <<' function just outputs the coords in // radians as I didn't want to drag in SLALIB. Instead, I use: // Astrometry::raToStr() & Astrometry::decToStr() in // echAstrometry/astrometry.hh e << "Supplied Object Position (app) RA = " << std::setprecision(6) << ra << "(" << ra_sign << std::setw(2) << hmsf[0] << ":" << hmsf[1] << ":" << hmsf[2] << "." << hmsf[3] << ") Dec = " << dec << "(" << dec_sign << dmsf[0] << ":" << dmsf[1] << ":" << dmsf[2] << "." << dmsf[3] << ")" << std::endl; sprintf(mesbuf, "Supplied Object Position (app) RA = %.6f(%c%.2d:%.2d:%.2d.%d) Dec = %.6f(%c%.2d:%.2d:%.2d.%d)", ra, ra_sign, hmsf[0], hmsf[1], hmsf[2], hmsf[3], dec, dec_sign, dmsf[0], dmsf[1], dmsf[2], dmsf[3]); } else e.dontThrow(); } /* *+ t e l M d l A p p 2 T a n * Function name: telMdl::App2Tan * Function: Determine Tangent plane coordinates from Apparent RA and Dec. * Description: Determines the Tangent coordinates on the focal plane, from the apparent RA and Dec of a source and the field center supplied to telMdl::CvtInit(). Accounts for refraction effects. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::App2Tan(ra, dec, xi, eta) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) ra (double) Apparent RA of source (radians). (>) dec (double) Apparent DEC of source (radians). (<) xi (double *) Tangent plane coordinate xi value (radians). (<) eta (double *) Tangent plane coordinate eta value (radians). * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ void telMdl::App2Tan ( double ra, double dec, double *pxi, double *peta) const { double ocra,ocdec; /* Observed ra,dec of field centre */ double ora,odec; /* Observed ra,dec of source */ double oaz,ozd,oha; double xi,eta; double pa; /* position angle of y axis */ int jstat; telMdl::ValidateApp2Tan(ra, dec); /* * Get the observed positions of the field center and source. Note the * different aoprms used - to allow for the telescope being pointed at a * different wavelength than we are observing at. * * Observed position of field centre */ slaAopqk(cenRa, cenDec, const_cast(cenAoprms), &oaz, &ozd, &oha, &ocdec, &ocra); /* * Observed position of source */ slaAopqk(ra, dec, const_cast(obsAoprms), &oaz, &ozd, &oha, &odec, &ora); /* * Project to tangent plane. Remember to catch errors. */ slaDs2tp(ora, odec, ocra, ocdec, &xi, &eta, &jstat); if (jstat != 0) THROW("Slalib returns error code %d from slaDs2tp",jstat); /* * Flip sign of xi for east negative */ xi = -xi; slaAopqk(ra, dec, const_cast(cenAoprms), &oaz, &ozd, &oha, &odec, &ora); /* * Rotate through ypa */ // XXX: The correct formula for rotation is: // x' = x . cos Ø - y . sin Ø // y' = x . sin Ø + y . cos Ø // SLALIB appears to use a different (non-conventional?) reference // frame: +y axis is _down_ (or rotation is clockwise, depending on // your paradigm.) // We need to be wary of this difference. pa = ypa; *pxi = xi * cos(pa) + eta * sin(pa); *peta = - xi * sin(pa) + eta * cos(pa); } /* * Class Private Function name: telMdl::App2Xy * Description: Determine the x and y coorindates on the field plate focal place from the apoparent RA and Dec of a soruce. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) ra (double) Apparent RA of source. (>) dec (double) Apparent DEC of source. (<) x (double *) X position of source on field plate (microns). (<) y (double *) Y position of source on field plate (microns). * History: 20-Feb-2003 - TJF - Original version */ void telMdl::App2Xy ( double ra, double dec, double *x, double *y) const { double xi,eta; /* tangent plane coordinates */ /* * Convert to tangent plane coorindates. */ App2Tan(ra, dec, &xi, &eta); /* * Apply Distortion convert to x/y. */ Tan2Xy(xi, eta, x, y); } /* * Class Private Function name: telMdl::Xy2Pos * Description: Convert the x/y coordinates of a object on the field plate, as output by telMdlApp2Xy to the coordinates required for input to the Echidna robot positioner system. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) x (double) Field plate x coordinate of object (microns). (>) y (double) Field plate y coordinate of object (microns). (>) *xp (double) Positioner X coordinate of object (microns). (>) *yp (double) Positioner Y coordinate of object (microns). * History: 20-Feb-2003 - TJF - Original version 19-Dec-2005 - SS - Comment. */ void telMdl::Xy2Pos ( double x, double y, double *xp, double *yp) const { double a,b,c,d,e,f; /* * Extract the transformation coefficent details. */ a = linPars.coeffs[0]; b = linPars.coeffs[1]; c = linPars.coeffs[2]; d = linPars.coeffs[3]; e = linPars.coeffs[4]; f = linPars.coeffs[5]; /* * Normalize the values. */ NormalizeLinear(linPars.rotation, 1, linPars.nonperp, &a, &b, &c, &d, &e, &f); /* * Apply the transformation. This is effectively doing a slaXy2xy(). */ *xp = a + b*x + c*y; *yp = d + e*x + f*y; } /* * Class Private Function name: telMdl::Xy2App * Description: Determine the apparent RA an Dec of a source from the x,y coorindate in the field plate and the field center. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) x (double) X position of source on field plate (microns). (>) y (double) Y position of source on field plate (microns). (<) ra (double *) Apparent RA of source. (<) dec (double *) Apparent DEC of source. * History: 20-Feb-2003 - TJF - Original version */ void telMdl::Xy2App ( double x, double y, double *ra, double *dec) const { double mra,mdec; /* Observed ra,dec of source */ double xi,eta; /* tangent plane coordinates */ double oaz; /* Observed Az */ double ozd; /* Observed Zd */ double oha; /* Observed HA */ double ocra,ocdec ; /* Observed field centre RA and Dec */ double pa; double xi2,eta2; /* * Observed position of field centre */ slaAopqk(cenRa,cenDec,const_cast(cenAoprms),&oaz,&ozd,&oha, &ocdec,&ocra); /* * Convert from focal plane to tanget plane, removing distortion. */ Xy2Tan(x, y, &xi, &eta); /* * Rotate through -ypa */ pa = -(ypa); xi2 = xi * cos(pa) + eta * sin(pa); eta2 = - xi * sin(pa) + eta * cos(pa); /* * Flip sign of xi */ xi2 = -xi2; /* * Deproject from tangent plane to give observed ra, dec */ slaDtp2s(xi2,eta2,ocra,ocdec,&mra,&mdec); /* * Observed to Apparent */ slaOapqk("R",mra,mdec,const_cast(obsAoprms),ra,dec); } /* * Class Private Function name: telMdl::Pos2Xy * Description: Convert the x,y coorindates of a star on the coorindate system used by the Echidna positioner task to the x,y coordinates on the field plate as used by Xy2App. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) xp (double) Positioenr X coordinate of object (microns). (>) yp (double) Positioenr Y coordinate of object (microns). (>) *x (double) Field plate x coordinate of object (microns). (>) *y (double) Field plate y coordinate of object (microns). * History: 20-Feb-2003 - TJF - Original version 23-Jan-2008 - TJF - Should check scale against 0, but instead check for a range around zero - since floating point compare to 0 is not certain to work. 25-Jan-2008 - TJF - Use Scott's FEQUAL macros to check floating point equality. */ void telMdl::Pos2Xy ( double xp, double yp, double *x, double *y) const { double a,b,c,d,e,f; /* * Extract inversion transformation coefficents from the model details. */ a = linPars.invCoeffs[0]; b = linPars.invCoeffs[1]; c = linPars.invCoeffs[2]; d = linPars.invCoeffs[3]; e = linPars.invCoeffs[4]; f = linPars.invCoeffs[5]; if (FEQUAL3(linPars.scale, 0)) THROW("Illegal scale. Must not be zero"); NormalizeLinear(-1.0*linPars.rotation, 1.0, -1.0*linPars.nonperp, &a, &b, &c, &d, &e, &f); /* * Apply the transformation. */ *x = a + b*xp + c*yp; *y = d + e*xp + f*yp; } /*^L *+ t e l M d l : : Q u i c k A p p t o O b s * Function name: telMdl::QuickAppToObs * Function: Perform a quick apparent to observed place conversion. * Description: Perform a quick apparent to observed place conversion using slaAopqk and parameters supplied to telMdl::CvtInit(). * Language: C++ * Namespace:echidna * Call: (void) = telMdl::QuickAppToObs(ra, dec, aob, zob, hob, dob, rob, ) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) ra (double) Apparent center RA (radians). (>) dec (double) Apparent Center Declination (radians). (<) aob (double *) Observed azimuth (radians, N=0, E = 90degrees) (<) zob (double *) Observed zenith distance (radians) (<) hob (double *) Observed Hour Angle (radians) (<) dob (double *) Observed Declination (radians) (<) rob (double *) Observed Right Ascension (radians) * Include files: echTelMdl.h *- * External functions used: * Support: Tony Farrell, AAO * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ void telMdl::QuickAppToObs ( double ra, double dec, double *aob, double *zob, double *hob, double *dob, double *rob) const { slaAopqk(ra,dec,const_cast(obsAoprms),aob,zob,hob,dob,rob); } /*^L *+ t e l M d l : : N o r m a l i z e * Function name: telMdl::Normalize * Function: Normalize our model. * Description: Apply extra rotation and non-perpendicularity effects to the linear model parameters and scale effects the distortion parameters. The extra parameters are then set to their null values (one for scale and zero for the other two). Note, rotation and non-perpendicularity are in degrees. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::Normalize() * Include files: echTelMdl.h *- * Support: Tony Farrell, AAO * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ void telMdl::Normalize () { /* * Normalize any rotation through the linear model. */ NormalizeLinear(linPars.rotation, 1, linPars.nonperp, &linPars.coeffs[0], &linPars.coeffs[1], &linPars.coeffs[2], &linPars.coeffs[3], &linPars.coeffs[4], &linPars.coeffs[5]); /* * Normalize any scale through the distortion model. */ NormalizeDist(linPars.scale, &distPars[0], &distPars[1], &distPars[2], &distPars[3], &distPars[4], &distPars[5]); /* * If sucessfull (no exception) - put these back to their default values. */ linPars.rotation = 0; linPars.scale = 1; linPars.nonperp = 0; } /* * Class Private Function name: telMdl::NormalizeLinear * Description: Apply rotation and scale effects to linear model parameters. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) rotation (double) A rotation in degrees. (>) scale (double) A scale (>) nonperp (double) Non-perpendicularity in degrees (!) a (double) Linear model parameter a (!) b (double) Linear model parameter b (!) c (double) Linear model parameter c (!) d (double) Linear model parameter d (!) e (double) Linear model parameter e (!) f (double) Linear model parameter f * History: 20-Feb-2003 - TJF - Original version 23-Jan-2008 - TJF - Should check scale against 0, but instead check for a range around zero - since floating point compare to 0 is not certain to work. 25-Jan-2008 - TJF - Use Scott's FEQUAL macros to check floating point equality. */ void telMdl::NormalizeLinear ( double rotation, double scale, double nonperp, double *a, double *b, double *c, double *d, double *e, double *f) { if (FEQUAL3(scale,0)) THROW("Illegal scale. Must not be zero"); else if (!FEQUAL9(scale,1)) // Checking down to 9 decimal places { *b *= scale; *c *= scale; *e *= scale; *f *= scale; } if (!FEQUAL9(rotation,0)) // Checking down to 9 decimal places { double cosine; double sine; double b2, c2, e2, f2; double an, bn, cn, dn, en, fn; cosine = cos(rotation/180.0*M_PI); sine = sin(rotation/180.0*M_PI); /* * Short hand for rotation coefficents */ b2 = cosine; c2 = sine; e2 = -1 * sine; f2 = cosine; /* * The original linear model is * x1 = a + b.x0 + c.y0 * y1 = d + e.x0 + f.y0 * * We want to apply this. * * x2 = b2.x1 + c2.y1 * y2 = e2.x1 + f2.y1 * * Which is equivalent to this. * * x2 = b2.(a + b.x0 + c.y0) + c2.(d + e.x0 + f.y0) * y2 = e2.(a + b.x0 + c.y0) + f2.(d + e.x0 + f.y0) * * And working out the new model, we get. * x2 = (b2.a + c2.d) + (b2.b + c2.e) x0 + (b2.c + c2.f) y0 * y2 = (e2.a + f2.d) + (e2.b + f2.e) x0 + (e2.c + f2.f) y0 * * Hence we can see our new parameters. */ an = b2*(*a) + c2*(*d); bn = b2*(*b) + c2*(*e); cn = b2*(*c) + c2*(*f); dn = e2*(*a) + f2*(*d); en = e2*(*b) + f2*(*e); fn = e2*(*c) + f2*(*f); *a = an; *b = bn; *c = cn; *d = dn; *e = en; *f = fn; } if (!FEQUAL9(nonperp,0)) { double cosine; double sine; double b2, c2, e2, f2; double an, bn, cn, dn, en, fn; double radians = nonperp/180*M_PI; cosine = cos(radians/2.0); sine = sin(radians/2.0); /* * * * A non-perpendicularity transformation is given by * * x' = cos(perp/2)*x + sin(perp/2)*y * y' = sin(perp/2)*x + cos(perp/2)*y * * Short hand for coefficents. */ b2 = cosine; c2 = sine; e2 = sine; f2 = cosine; /* * The original linear model is * x1 = a + b.x0 + c.y0 * y1 = d + e.x0 + f.y0 * * We want to apply this. * * x2 = b2.x1 + c2.y1 * y2 = e2.x1 + f2.y1 * * Which is equivalent to this. * * x2 = b2.(a + b.x0 + c.y0) + c2.(d + e.x0 + f.y0) * y2 = e2.(a + b.x0 + c.y0) + f2.(d + e.x0 + f.y0) * * And working out the new model, we get. * x2 = (b2.a + c2.d) + (b2.b + c2.e) x0 + (b2.c + c2.f) y0 * y2 = (e2.a + f2.d) + (e2.b + f2.e) x0 + (e2.c + f2.f) y0 * * Hence we can see our new parameters. */ an = b2*(*a) + c2*(*d); bn = b2*(*b) + c2*(*e); cn = b2*(*c) + c2*(*f); dn = e2*(*a) + f2*(*d); en = e2*(*b) + f2*(*e); fn = e2*(*c) + f2*(*f); *a = an; *b = bn; *c = cn; *d = dn; *e = en; *f = fn; } } /* * Class Private Function name: telMdl::NormalizeDist. * Description: Apply a scale effect to the distortion model parameters. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) scale (double) A scale (!) d0 (double) Distortion model parameter 0 (!) d1 (double) Distortion model parameter 1 (!) d2 (double) Distortion model parameter 2 (!) d3 (double) Distortion model parameter 3 (!) d4 (double) Distortion model parameter 4 (!) d5 (double) Distortion model parameter 5 * History: 20-Feb-2003 - TJF - Original version 23-Jan-2008 - TJF - Should check scale against 0, but instead check for a range around zero - since floating point compare to 0 is not certain to work. 25-Jan-2008 - TJF - Use Scott's FEQUAL macros to check floating point equality. */ void telMdl::NormalizeDist ( double scale, double *d0, double *d1, double *d2, double *d3, double *d4, double *d5) { if (FEQUAL3(scale, 0)) THROW("Illegal scale. Must not be zero"); else if (!FEQUAL6(scale,1)) { *d0 *= scale; *d1 *= scale; *d2 *= scale; *d3 *= scale; *d4 *= scale; *d5 *= scale; } } /* * Class Private Function name: telMdl::GetFocalLength * Description: Determine the focal length and pin cushion distortion effect parameter. These allow us to do an approximate distortion conversion. * History: 20-Feb-2003 - TJF - Original version 25-Jan-2008 - TJF - Due to changes in the distortion calculation, this changes. Comments updated and improved. */ void telMdl::GetFocalLength( double *focalLength, double *pincushionCoeff) const { const double d0 = distPars[0]; const double d1 = distPars[1]; /* * Using only first two parameters of the distortion model, we have * * R = d0 * th + d1 * th^3 * * * th is the radius to the object from the field center in degrees, * l is the observation wavelength in nm, * R is the field plate radius meters. * d0 is m/degree * d1 is (m/degree)^3 * * The focal length in meters FL = d0 * 180/pi. * * Let Ra = R/FL * * Ra = (d0 * th + d1 * th^3)/FL * * Let r = (th/180*pi), hence th = r*180/pi * Thus: * * Ra = ((d0 * r * 180/pi) + (d1 * r^3 * 180^3 / pi^3))/FL * = ( r ( d0 * 180/pi + d1 * 180^3 / pi^3 * r^2))/FL * = r ( 1 + d1/d0 * 180^2/pi^2 * r^2) * * And this is of the same format at the pincushion effect formula - see * for example, the slaPcd() documentation, the formula is * p = r(1 + c r^2) * * * As a result, our pincushion coffecient is * d1/d0 * (180/pi)^2 * * Remember that we want the focal length in microns rather then m */ *focalLength = d0 * 180 / M_PI * 1000000.0 * linPars.scale; *pincushionCoeff = d1 * pow((180.0/M_PI), 2) / d0; } /* * Internal Function name: t e l M d l : : a p p l y C o r r e c t o r D i s t o r t i o n * Description: Calculate the position at which point $r is distorted to by the corrector. * History: 10-Jun-2004 - SS - Initial version; extracted from telMdl::Tan2Xy(). 10-Oct-2006 - SS - Use $TEL_MDL_FOV_WARN for threshold value. 08-Feb-2007 - SS - Created $USE_OLD_DIST_MDL for debugging telescope model with OzPos data. 22-Jan-2008 - TJF - Historical note on previous version - The new formula was correct against the spreadsheet, I (TJF) has misread the spreadsheet. Also - you need to use the values from the spreadsheet, getting you values in meters, rather then trying to multiple them, otherwise the values won't fit. The code used in the previous formula corresponds to spreadsheet sheet named "coded-old-version". 22-Jan-2008 - TJF - Revert to original code as spreadsheet did not need to calculate the value (y=th*d0) and use it in the later terms. This would must make it harder to calibrate the various terms independelty and does not change the results. Use parameters and units directly from the spreadsheet. Any unit conversion might cause issues. */ double telMdl::applyCorrectorDistortion (const double th) const { double l = obsLambda * 1000.0; // Observation wavelength (nm) double d0, d1, d2, d3, d4, d5; // Distortion model params to be normalised static const double warnDist = getFloatEV("TEL_MDL_FOV_WARN", 0.60) / 2.0; if (th > warnDist) cerr << "Warning: calculating distortion " << th << "° from field centre" << endl; // Fetch the distortion parameters. d0 = distPars[0]; d1 = distPars[1]; d2 = distPars[2]; d3 = distPars[3]; d4 = distPars[4]; d5 = distPars[5]; // Normalize any scale into the distortion parameters. NormalizeDist(linPars.scale, &d0, &d1, &d2, &d3, &d4, &d5); /* * This formula was determined by Peter Gillingham and Tony Farrell. * See the spreadsheet distortion_calc2.xls, page "as-per-new-code". PG * developed the "pg_orig" page and TJF added the other pages to develop this * formula. * * The best fit calcuated by the Excel solver is given by . * * R = 0.29195 th + 0.0445 th^3 + 0.0364 th^5 + (-0.197) th/l * * This gives a max error of 2 microns and an RMS of 1.5 microns. * * th is the radius to the object from the field center in degrees, * l is the observation wavelength in nm, * R is the field plate radius meters. * * */ double R = (d0 * th ) + (d1 * pow(th,3) ) + (d2 * pow(th,5) ) + (d3 * th/l ) + (d4 * pow(th/l,2)); R *= 1000000.0; // Convert to micrometers #ifdef TESTCODE if (bDebugTM) { cout << "applyCorrectorDistortion(th=" << th << "°) = " << R << "µm" << endl; cout << "d0=" << d0 << " d1=" << d1 << " d2=" << d2 << " d3=" << d3 << " d4=" << d4 << " d5=" << d5 << " l=" << l << endl; } #endif return R; } /* * Class Private Function name: telMdl::Tan2XY * Description: Convert tangent plane coordindates to x/y taking account of distortion. * History: 20-Feb-2003 - TJF - Original version 04-Mar-2003 - TJF - Peter's first pass as distortion formula modifed by TJF used. 29-Aug-2003 - TJF - Modify Flexure model as per latest info. 12-Nov-2003 - SS - bugfixes. 10-Jun-2004 - SS - remove hardcoded constants added in last revision; extract code to create applyCorrectorDistortion(). 24-Jan-2006 - SS - Fix XY alignment of flexure compensation. 23-Jan-2008 - TJF - Should check scale against 0, but instead check for a range around zero - since floating point compare to 0 is not certain to work. Similarly for r. 24-Jan-2008 - TJF - Index flexure parameters using marcros instead of numbers. Get the distortion center and apply it to the calculations. 25-Jan-2008 - TJF - Use Scott's FEQUAL macros to check floating point equality. 27-Jan-2008 - TJF - Don't remove xZeroOffset/yZeroOffset since it as in effect already removed. Apply flexure compensation when x/y is zero. 28-Jan-2008 - TJF - flexure rotation calcuation used Y instead of X at one point. Also - use rotatorSetting_ for the rotator value rather then ypa, which is just the instrument offset. */ void telMdl::Tan2Xy ( const double xi, const double eta, double *x, double *y) const { /* * Fetch the offset of the distortion. */ double xiZeroOffset; double etaZeroOffset; double xZeroOffset; double yZeroOffset; GetDistCenter(&xiZeroOffset, &etaZeroOffset, &xZeroOffset, &yZeroOffset); /* * Determine the radius to the object, accounting for the * distortion centre offset. */ double r = hypot(xi - xiZeroOffset, eta - etaZeroOffset); // Tangent plane radius angle (radian.) if (FEQUAL3(linPars.scale, 0)) THROW("Illegal scale. Must not be zero"); // Avoid divide by zero errors. if (r > 0.00000001) { const double R = applyCorrectorDistortion(r / M_PI * 180.0); // R in (micron) // Get x and y components. *x = xi * R / r; // R/r represents (micron/radian) *y = eta * R / r; } else { *x = 0; *y = 0; } #ifdef TESTCODE // Determine the pincushion coefficent for slaPcd() and // the focal length determined by the distortion parameters. double focalLength; double pincushionCoeff; GetFocalLength(&focalLength, &pincushionCoeff); // Apply the distortion using slaPcd. Here we are just confirming // the new formula works ok and checking for degeneration. double xi1 = xi; double eta1 = eta; // Apply Radial Distortion slaPcd(pincushionCoeff, &xi1, &eta1); // Multiply by focal length to get x and y in µm. double pcdX = xi1 * focalLength; double pcdY = eta1 * focalLength; #ifdef TESTCODE if (bDebugTM) { cout << "*x=" << *x << endl; cout << "*y=" << *y << endl; cout << "xi=" << xi << endl; cout << "eta=" << eta << endl; cout << "pcdX=" << pcdX << endl; cout << "pcdY=" << pcdY << endl; cout << "UNpcdX=" << (xi * focalLength) << endl; cout << "UNpcdY=" << (eta * focalLength) << endl; cout << "coeff=" << pincushionCoeff << endl; } #endif #warning "Change required here - check against slaPCD result." // This should be change to an algorithm which checks the // slaPcd results and if that shows the result outside the // field of view, then we use those results rather then // the original distortion formula result. Alternatively, // try to work out where the distortion formula degenerates. #define MAX_DEVIATION 8000 /* Maximum deviation of proper formula from slaPcd formula within field of view. about 8mm edge */ /* * If the difference between the Radius worked out * this way and worked out using the proper formula * is more then MAX_DEVIATION, then it is possible * the proper forumula worked out by slaPcd. * * It is likely that th > 0.25 degrees in this case, * outside the ozpoz field of view. */ double deviation = fabs(hypot(pcdX, pcdY) - hypot(*x, *y)); if (deviation > MAX_DEVIATION) { fprintf(stderr,"Warning - object position with model deviated from PCD position\n"); fprintf(stderr," - Object deviation %6.1fµm (max dev. %dµm)\n", deviation, MAX_DEVIATION); #if 0 fprintf(stderr," - Using position generated by slaPcd() instead of model\n"); *x = pcdX; *y = pcdY; #endif } #endif //Ensure that we don't allow really high values. if (*x > X_Y_MAX || *y > X_Y_MAX) THROW("FP coords (%d, %d) are too large.", x, y); /* * Flexure compensation. * * Taken from an e-mail from Masayuki Akiyama sent on 26-aug-2003 to Scott * Smedley). * * TODO: Why is it necessary to have coefficients in mm? * * The coefficent is appropiate for shift in mm. - so must multiply * by 1000 to get µm. * * telMdlPAR_FLEX_FIRST is an index into the telMdlPARS structure. We need to * index into the distPars variable. Element 0 of this is equivalent * to telMdlPAR_NW_DIST_FIRST in telMdlPARS, so we subtract this to get * out index. */ double fc1 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST]; double fc2 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST+1]; double fc3 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST+2]; double fc4 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST+3]; double el = M_PI/2 - zenithDistance; // Elevation double flexShiftX = (fc1 * sin(el) + fc2 * cos(el)) * 1000.0; double flexShiftY = (fc3 * sin(el) + fc4 * cos(el)) * 1000.0; #ifdef TESTCODE if (bDebugTM) { cout << "deviation=" << deviation << endl; cout << "flexShiftX=" << flexShiftX << " flexShiftY=" << flexShiftY << endl; } #endif /* * $flexShiftY indicates the amount of offset in the pseudo-vertical * direction. The instrument Y axis is likely different from this, * so we need to rotate appropriately. * TODO: Should we be using the SLALIB-type rotation which (effectively) * rotates in the opposite direction? (See comments in App2Tan().) * In fact, which direction should we rotated in? * */ *x += flexShiftX * cos(rotatorSetting_) - flexShiftY * sin(rotatorSetting_); *y += flexShiftX * cos(rotatorSetting_) + flexShiftY * sin(rotatorSetting_); } /* * Class Private Function name: telMdl::X Y 2 T a n * Description: Convert x/y to tangent plane coordinates taking account of distortion. Note - this conversion is approximate as we don't have a real function for doing this. This should be ok was this is only used for determining local guiding offset and should be locally correct to the required precision. * History: 20-Feb-2003 - TJF - Original version 03-Mar-2003 - TJF - Initial flexure correction added. 29-Aug-2003 - TJF - Modify Flexure model as per latest info. 24-Jan-2006 - SS - Fix XY alignment of flexure compensation. 23-Jan-2008 - TJF - Should check scale against 0, but instead check for a range around zero - since floating point compare to 0 is not certain to work. 24-Jan-2008 - TJF - Index flexure parameters using marcros instead of numbers. Get the distortion center and apply it to the calculations. 25-Jan-2008 - TJF- In making the last change, I had created some new variables (xi1, eta1) and then not initialised them. Fixed by version 2.1 25-Jan-2008 - TJF - Use Scott's FEQUAL macros to check floating point equality. 27-Jan-2008 - TJF - Remove distoration center stuff and add comment suggesting this does not work anyway and we are saved by the iterative solution. 28-Jan-2008 - TJF - flexure rotation calcuation used Y instead of X at one point. Also - use rotatorSetting_ for the rotator value rather then ypa, which is just the instrument offset. */ void telMdl::Xy2Tan ( double x, double y, double *xi, double *eta) const { double focalLength; double pincushionCoeff; if (FEQUAL3(linPars.scale, 0)) THROW("Illegal scale. Must not be zero"); /* * Remove flexure compensation. (Remember, coefficent in mm, we work * in µm). */ double fc1 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST]; double fc2 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST+1]; double fc3 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST+2]; double fc4 = distPars[telMdlPAR_FLEX_FIRST-telMdlPAR_NW_DIST_FIRST+3]; double el = M_PI/2 - zenithDistance; // Elevation // TODO: Why do flexure compensation coefficient have to be in mm? double flexShiftX = (fc1 * sin(el) + fc2 * cos(el)) * 1000.0; double flexShiftY = (fc3 * sin(el) + fc4 * cos(el)) * 1000.0; x -= flexShiftX * cos(rotatorSetting_) - flexShiftY * sin(rotatorSetting_); y -= flexShiftX * cos(rotatorSetting_) + flexShiftY * sin(rotatorSetting_); if (x > X_Y_MAX || y > X_Y_MAX) THROW("FP coords (%d, %d) are too large.", x, y); // // Note - this is suspect - see tdfxy.c/TdfDistXyInv(), but // may only be being used in the iterative code - so may // not be huring us. // // // Determine the focal length and pincushion coefficent from the // telescope model parameters. GetFocalLength(&focalLength, &pincushionCoeff); *xi = x / focalLength; *eta = y / focalLength; // Remove Radial Distortion (approximate only) slaUnpcd(pincushionCoeff, xi, eta); } /* *+ t e l M d l : : M i c r o n s 2 A r c s e c * Function name: telMdl::Microns2Arcsec * Function: Convert µm on the focal plane to arc-seconds on the sky. * Description: Uses the model to convert the specified number of microns on the focal plane to approximate arc-seconds on the sky (most accurate near the plate center). * Language: C++ * Namespace:echidna * Call: (double) = telMdl::Microns2Arcsec(microns) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) microns (double) The microns to convert * Returned value: The equivalent arc-seconds at the plate centre. * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ #define SEC2RAD(_sec) ((_sec)/(180.0*60.0*60.0)*M_PI) extern double telMdl::Microns2Arcsec( const double microns) const { double x, y; double rx, ry; double arcsec2microns; /* * Convert 1 arc second to X and Y. */ Tan2Xy(0, SEC2RAD(1.0), &x, &y); /* * Now convert to robot coordindate system. */ Xy2Pos(x, y, &rx, &ry); /* * Now we can work out how many microns in an arcsecond (at least * clost to the plate centre) */ arcsec2microns = sqrt(rx*rx + ry*ry); return (microns/arcsec2microns); } /* *+ t e l M d l : : A p p 2 P o s * Function name: telMdl::App2Pos * Function: Determine Echidna coordinates from Apparent RA and Dec. * Description: Determines the Echdina robot coordinates from the apparent RA and Dec of a source and field center supplied to telMdl::CvtInit(). Performs a full conversion accounting for refraction and distortion. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::App2Pos(ra, dec, x, y) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) ra (double) Apparent RA of source (radians). (>) dec (double) Apparent DEC of source (radians). (<) x (double *) Echidna FPI coordinates x value (microns). (<) y (double *) Echidna FPI coordinates y value (microns). * Include files: echTelMdl.h *- FUNCTION DEFINED IN HEADER FILE - HEADER FOR DOC GENERATION ONLY. * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ /* *+ t e l M d l : : P o s 2 A p p * Function name: telMdl::Pos2App * Function: Determine approximate Apparent RA and Dec from Echidna x/y. * Description: Determines the approximate Apparent RA and Dec from echidna x/y using values supplied to telMdl::CvtInit(). The conversion is only approximate, as the forward distortion model does not have a known inverse. Pincushion distortion is assumed and slaUnpcd() used to perform the removal of the distortion. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::Pos2App(x, y, ra, dec) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) x (double) Echidna FPI coordinates - x value (microns). (>) y (double) Echidna FPI coordinates - y value (microns). (<) ra (double *) Apparent RA of source (radians). (<) dec (double *) Apparent DEC of source (radians). * Include files: echTelMdl.h *- FUNCTION DEFINED IN HEADER FILE - HEADER FOR DOC GENERATION ONLY. * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ /* *+ t e l M d l : : P o s 2 A p p I t e r * Function name: telMdl::Pos2AppIter * Function: Determine _precise_ Apparent RA and Dec coords for a given set of field plate coords. * Description: The telMdl::Pos2App() member function _approximates_ the apparent RA/Dec for a given set of field plate coords. This member function aims to calculate a more accurate value for the apparent RA/Dec. This is accomplished by checking the return result of Pos2App() against its inverse - App2Pos(). If the field plate coords disagree by more than $tolerance, this routine will adjust the RA/Dec coords & try again. Up to $maxIters iterations may occur before this function returns. * Language: C++ * Namespace:echidna * Call: (void) = telMdl::Pos2AppIter(x, y, ra, dec) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) x (double) Echidna FPI coordinates - x value (microns). (>) y (double) Echidna FPI coordinates - y value (microns). (<) ra (double *) Apparent RA of source (radians). (<) dec (double *) Apparent DEC of source (radians). (>) tolerance (double *) accuracy of $ra, $dec coords on field plate. (>) maxIters (double *) maximum number of iterations. * Include files: echTelMdl.h * See Also: telMdlMean::Pos2App() member function. * Support: Scott Smedley, AAO *- * History: 17-Jun-2004 - SS - Initial version. */ void telMdl::Pos2AppIter ( const double x, const double y, double *ra, double *dec, const double tolerance, const int maxIters) const { Pos2App(x, y, ra, dec); const double initialRA = *ra, initialDec = *dec; for (int i = 1; i <= maxIters; ++i) { double nx, ny; App2Pos(*ra, *dec, &nx, &ny); const double deltaX = nx - x; const double deltaY = ny - y; const double delta = hypot(deltaX, deltaY); if (delta <= tolerance) break; double nRA, nDec, deltaRA, deltaDec; Pos2App(x + deltaX, y + deltaY, &nRA, &nDec); deltaRA = initialRA - nRA; deltaDec = initialDec - nDec; *ra += deltaRA; *dec += deltaDec; if (i == maxIters) cerr << "Warning: reached $maxIters in Pos2AppIter()" << endl; } } /* *+ t e l M d l : : E q u * Function name: telMdl::DeltaEqu * Function: Determine delta x-y/ra-dec for a given Apparent RA/Dec within the field. * Description: This method determines the change in RA & Dec (apparent) indicated by a given change in x &y. The function returns factors such that DeltaRa = DeltaX * drdx + DeltaY * drdy DeltaDec = DeltaY * dddx + DeltaY * dddy This conversion will be local to the specified RA/Dec position (it can't be used in other locations due to field distortion) * Call: * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) appRa (const double) The RA at which we determine the conversion. Radians. (>) appDec (const double) The Dec at which we determine the conversion. (<) drdx (double *) Change in RA for a change in X (<) drdy (double *) Change in RA for a change in Y (<) dddx (double *) Change in Dec for a change in X (<) dddy (double *) Change in Dec for a change in Y * Returned value: * See Also: telMdlMean::DeltaEqu(), telMdl::DeltaAltAz(). * Support: Tony Farrell, AAO *- * History: 23-Jan-2008 - THF - Initial version. */ const double OFFSET = DAS2R*5; void telMdl::DeltaEqu( const double appRa, const double appDec, double * const drdx, /* Change in RA for a change in x */ double * const drdy, /* Change in RA for a change in y */ double * const dddx, /* Change in Dec for a change in x */ double * const dddy /* Change in Dec for a change in y */ ) const { double baseX; double baseY; // Get the position of the object on the focal plane) App2Pos(appRa, appDec, &baseX, &baseY); // Move RA by the specified offset and calculate new X and Y // QUESTION - should we use cos(dec) here? double raOffX, raOffY; App2Pos(appRa+OFFSET, appDec, &raOffX, &raOffY); // Move Dec by the specified offset and caculate new X and Y double decOffX, decOffY; App2Pos(appRa, appDec+OFFSET, &decOffX, &decOffY); // Caclualte the results. Units are arc-seconds/micron *drdx = OFFSET/(raOffX - baseX); *drdy = OFFSET/(raOffY - baseY); *dddx = OFFSET/(decOffX - baseX); *dddy = OFFSET/(decOffY - baseY); } /* *+ t e l M d l : : E q u * Function name: telMdl::DeltaAltAz * Function: Determine delta x-y/az-el for a given Apparent RA/Dec within the field. * Description: This method determines the change in Azimuth and Elevation indicated by a given change in x &y. The function returns factors such that DeltaAz = DeltaX * dadx + DeltaY * dady DeltaEl = DeltaY * dedx + DeltaY * dedy This conversion will be local to the specified mean RA/Dec position (it can't be used in other locations due to field distortion). Note - we do mean the specified RA/Dec position, not an Alt/Az position. * Call: * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) appRa (const double) The RA at which we determine the conversion. Radians. (>) appDec (const double) The Dec at which we determine the conversion. (<) dadx (double *) Change in Azimuth for a change in X (<) dady (double *) Change in Azimuth for a change in Y (<) dedx (double *) Change in RA for a change in X (<) dedy (double *) Change in RA for a change in Y * Returned value: * See Also: telMdlMean::DeltaAltAz(), telMdl::DeltaEqu(). * Support: Tony Farrell, AAO *- * History: 23-Jan-2008 - THF - Initial version. */ void telMdl::DeltaAltAz( const double appRa, const double appDec, double * const dadx, /* Change in Azimuth for a change in x */ double * const dady, /* Change in Azimuth for a change in y */ double * const dedx, /* Change in Elevation for a change in x */ double * const dedy /* Change in Elevation for a change in y */ ) const { double baseX; double baseY; double az; /* Azimuth at this position */ double el; /* Elevation at this position */ Eq2AltAz(appRa, appDec, &az, &el); // Get the position of the object on the focal plane) //App2Pos(appRa, appDec, &baseXchk, &baseYchk); AltAz2Pos(az, el, &baseX, &baseY); // Move Azimuth by the specified offset and calculate new X and Y double azOffX, azOffY; AltAz2Pos(az+OFFSET, el, &azOffX, &azOffY); // Move Elevation by the specified offset and caculate new X and Y double elOffX, elOffY; App2Pos(az, el+OFFSET, &elOffX, &elOffY); // Caclualte the results. Units are arc-seconds/micron *dadx = OFFSET/(azOffX - baseX); *dady = OFFSET/(azOffY - baseY); *dedx = OFFSET/(elOffX - baseX); *dedy = OFFSET/(elOffY - baseY); } /* * Class Private Function name: telMdl::AltAz2Xy * Description: Determine the x and y coorindates on the field plate focal place from the Azimuth and elevation of a source. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) az (double) Azimuth of source. (>) el (double) Elevation DEC of source. (<) x (double *) X position of source on field plate (microns). (<) y (double *) Y position of source on field plate (microns). * History: 23-Jan-2008 - TJF - Original version */ void telMdl::AltAz2Xy ( double az, double el, double *x, double *y) const { double xi; double eta; /* * Convert to tangent plane coorindates. */ AltAz2Tan(az, el, &xi, &eta); /* * Apply Distortion convert to x/y. */ Tan2Xy(xi, eta, x, y); } /* *+ t e l M d l : : E q 2 A l t A z * Function name: telMdl::Eq2AltAz * Function: Determine Azimuth and Elevation of an object in the field given the equatorial position. * Description: Uses slaAopqk() and the current apparent to obseved conversion parameters. * Call: * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) appRa (const double) The RA Radians. (>) appDec (const double) The Dec Radians. (<) az (double *) The azimuth (radians) (<) el (double *) The elevation (radians) * Returned value: * See Also: * Support: Tony Farrell, AAO *- * History: 23-Jan-2008 - THF - Initial version. */ void telMdl::Eq2AltAz( double ra, double dec, double *az, double *el) const { double ozd; double oha; double ocdec; double ocra; /* * slaAopqk does the for us - but note that it returns * zd intead of elevation. */ slaAopqk(ra, dec, const_cast(obsAoprms), az, &ozd, &oha, &ocdec, &ocra); /* * Work out the elevation. */ *el = DPIBY2 - ozd; } /* *+ t e l M d l : : A l t A z 2 T a n * Function name: telMdl::AltAz2Tan * Function: Determine Tangent plane coordinates from Azimuth and Elevation. * Description: Determines the Tangent coordinates on the focal plane, from the Azimuth and Elevation of a source and the field center supplied to telMdl::CvtInit(). Accounts for refraction effects. Annonlying - at the moment this converts back to apparent Ra/Dec and then uses App2Tan() to do the work - to avoid repeating or breaking up App2Tan(). * Language: C++ * Namespace:echidna * Call: (void) = telMdl::App2Tan(ra, dec, xi, eta) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) ra (double) Apparent RA of source (radians). (>) dec (double) Apparent DEC of source (radians). (<) xi (double *) Tangent plane coordinate xi value (radians). (<) eta (double *) Tangent plane coordinate eta value (radians). * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ void telMdl::AltAz2Tan( double az, double el, double *xi, double *eta ) const { double ra, dec; /* Convert from az/el to apparent ra/dec */ slaOapqk ("A", az, el, const_cast(obsAoprms), &ra, &dec); /* Now convert to the tangent plane */ App2Tan(ra, dec, xi, eta); } /* * Class Private Function name: telMdl::GetDistCenter * Description: Determine the center of the distortion on the tanget plane. The distortion model parameters to will contain 1. The distortion offset in microns in x and y. 2. The RA/Dec position of the field at which this was true. 3. The Time at which this field was observed. 4. Any Position offset which had been applied to the rotator. From 2 and 3 above we can calculate the parallactic angle associated with the obsevation. From this, we remove the rotator position offset (4) to get a "rotator setting" at which the x/y (1) value was determined. From the current field information, we determine the parallatic angle of the current observation. We add in any rotator position offset to get a "rotator setting" for the current field. We rotate the x and y values (1) by the difference between the two rotator settings, to get the offset in x and y at the current rotator setting. We then divide this by the rough plate scale to get the xi and eta values of interest. * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (<) xi (double *) center in xi. (<) eta (double *) center in eta (<) x (double *) center in x (microns) (<) y (double *) center in y (microns) * History: 24-Jan-2008 - TJF - Original version */ typedef struct CenterOffCalDataType { double x0; double y0; double zeroRA; double zeroDEC; double zeroMJD; double zeroPA; }; void telMdl::GetDistCenter( double * const xi, double * const eta, double * const x, double * const y) const { *xi = 0; *eta = 0; *x = 0; *y = 0; // return; /* * The parameters last time we were run. We need to recalculate if they change, * but don't want to recalculate if they have not changed, so store the values * in a static structure. */ static CenterOffCalDataType savedData={0,0,0,0,0,0}; /* * The values we want to return. */ static double zeroPntXi = 0; static double zeroPntEta = 0; static double zeroPntX = 0; static double zeroPntY =0; /* * telMdlPAR_ZOFF_X etc are indices into the telMdlPARS structure. We need to * index into the distPars variable. Element 0 of this is equivalent * to telMdlPAR_NW_DIST_FIRST in telMdlPARS, so we subtract this to get * out index. */ CenterOffCalDataType curData; /* For storing the parameters now */ curData.x0 = distPars[telMdlPAR_ZOFF_X-telMdlPAR_NW_DIST_FIRST]; curData.y0 = distPars[telMdlPAR_ZOFF_Y-telMdlPAR_NW_DIST_FIRST]; curData.zeroRA = distPars[telMdlPAR_ZOFF_RA-telMdlPAR_NW_DIST_FIRST]; curData.zeroDEC = distPars[telMdlPAR_ZOFF_DEC-telMdlPAR_NW_DIST_FIRST]; curData.zeroMJD = distPars[telMdlPAR_ZOFF_MJD-telMdlPAR_NW_DIST_FIRST]; curData.zeroPA = distPars[telMdlPAR_ZOFF_PA-telMdlPAR_NW_DIST_FIRST]; /* * If savedData and curData are not identical, then savedData is * invalid and we must calcuate "savedRotation". We do it this way because * we are called many times and would not normally see any of these * parameters changed. */ if (memcmp((void *)&savedData, (void *)&curData, sizeof(savedData)) != 0) { /* * The temp, pressure and humidity dont significantly affect the * calculation of the parallactic angle - so use default values. */ const double TEMP=273; const double PRES=1013; const double HUMID=0.1; /* * Save the new parametes */ savedData = curData; /* * Convert the mean position to an apparent position. */ double appRa, appDec; slaMap(curData.zeroRA, curData.zeroDEC, 0, 0, 0, 0, 2000, curData.zeroMJD, &appRa, &appDec); /* * Now convert to an observed position. */ double oaz,ozd,oha, odec, ora; slaAop(appRa, appDec, curData.zeroMJD, dut,-longit_,lat_,height_,0.0,0.0, TEMP, PRES, HUMID, cenLambda, TRP_LSP_RATE, &oaz,&ozd,&oha,&odec,&ora); /* * Determine the rotation at this position - which is the parallactic angle (since * this is a prime focus instrument) */ double parallactic = slaPa(oha, odec, lat_); /* * The actual rotator position at this point includes the rotator position * angle */ double zeroPntRotation = parallactic - curData.zeroPA; /* * Rotate the distortion zero point by the difference between the * zero point rotator setting and our rotator setting. */ fprintf(stdout," rotatorSetting=%f, zeroPntRotation=%f \n", rotatorSetting_*180.0/M_PI, zeroPntRotation*180.0/M_PI); double rotatorDiff = rotatorSetting_ - zeroPntRotation; zeroPntX = curData.x0 * cos(rotatorDiff) - curData.y0 * sin(rotatorDiff); zeroPntY = curData.x0 * sin(rotatorDiff) + curData.y0 * cos(rotatorDiff); fprintf(stdout," Current Distortion Center FPIX = %f FPIY = %f \n", zeroPntX, zeroPntY ); /* * Scale to the tangent plane. The first coorindate in the distoriton * matrix is close enough, but its units are meters/degree. We want microns/radian */ double scale = distPars[0] * 1000000.0 * 180.0 / M_PI; // Convert to microns/radians. zeroPntXi = zeroPntX/scale; // radiauns zeroPntEta = zeroPntY/scale; // radiauns } /* * Return the values. */ *xi = zeroPntXi; *eta = zeroPntEta; *x = zeroPntX; *y = zeroPntY; } /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * call telMdlMean methods follow */ /* *+ t e l M d l M e a n : : A p p 2 P o s * Function name: telMdlMean::Mean2Pos * Function: Determine echidna coordinates from Mean RA and Dec. * Description: Determines the Echdina robot coordinates from the apparent RA and Dec of a source and field center supplied to telMdl::CvtInit(). * Language: C * Namespace:echidna * Call: (void) = telMdlMean::Mean2Pos(ra, dec, pmPars, x, y, appRa=0, appDec=0) const * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) ra (double) Mean RA of source (radians). (>) dec (double) Mean DEC of source (radians). (>) pmPars (const telMdlPM_PARS &) Proper motion parameters for source. (<) x (double *) Echidna FPI x coordinate (microns). (<) y (double *) Echidna FPI x coordinate (microns). (<) appRa (double *) If non-null, the apparent RA is written here (radians). (<) appDec (double *) If non-null, the apparent Dec is written here (radians). * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ void telMdlMean::Mean2Pos( const double meanRa, const double meanDec, const telMdlPM_PARS & pmPars, double * const robotX, double * const robotY, double * const appRa, double * const appDec) const { double ra; double dec; /* * Get apparent RA and Dec values, accounting for proper motion. */ slaMapqk (meanRa, meanDec, pmPars.ra, pmPars.dec, pmPars.parallax, pmPars.radialVel, const_cast(gmap), &ra, &dec); /* * If these are wanted, return them. */ if (appRa) *appRa = ra; if (appDec) *appDec = dec; /* * Translate sky-to-field positioner XY [microns] */ App2Pos (ra, dec, robotX, robotY); } /* *+ t e l M d l M e a n: : C v t I n i t * Function name: telMdlMean::CvtInit * Function: Initialize for conversions at a specified time, place etc. * Description: This methods initializes the class for conversions at particular time, wavelength and telescope model details and apparent place. This allows a sequence of transformation for one field to be done efficiently. The telMdlMean version of this function is handles mean telescope positions efficiently in conjunction with the Mean2Xy function. * Language: C++ * Namespace:echidna * Call: (void) = telMdlMean::CvtInit(mjd, dut, temp, pres, humid, cenWave, obsWave rotation, flexCorPars, ra, dec, equinox=2000) * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) mjd (double) UTC date and time as a modified Julian date. (>) dut (double) Delta UT (UT1 - UTC) seconds. (>) temp (double) Atmospheric temperature (K). (>) press (double) Atmospheric pressure (mB). (>) humid (double) Atmospheric humidity (0 to 1.0). (>) cenWave (double) Wavelength telescope is pointed with (microns) (>) obsWave (double) Wavelength of observation (micons). (>) rotation (double) Rotator position angle. (>) flexCorPars (const telMdlFlexCorPars &) flex correction parameters. (>) ra (double) Mean RA of observation. (>) dec (double) Mean Dec of observation. (>) equinox (double) Equniox of mean position. * Include files: echTelMdl.h *- * History: 12-Feb-2003 - TJF - Original Version, based on OzPoz version. */ void telMdlMean::CvtInit ( double mjd, double temp, double press, double humid, double cenWave, double obsWave, double rotation, const telMdlFlexCorPars &flexCorPars, double meanRa, double meanDec, double equinox) { double appRa; double appDec; // Set up for calling slaMapqk slaMappa(equinox, mjd, gmap); // Convert mean field centre RA/Dec to apparent/ (radians). slaMapqkz(meanRa, meanDec, gmap, &appRa, &appDec); telMdl::CvtInit(mjd, temp, press, humid, cenWave, obsWave, rotation, flexCorPars, appRa, appDec); } /* *+ t e l M d l M e a n : : P o s 2 M e a n * Function name: telMdlMean::Pos2Mean * Function: Convert field plate coordinates into mean RA/Dec coords. * Description: Note: this function returns approximate RA/Dec coords only. This is due to the inverse distortion model calculations being an approximation. Use telMdlMean::Pos2MeanIter() for a more accurate RA/Dec value. * Call: * Parameters: (">" input, "!" modified, "W" workspace, "<" output) > x (const double) x field plate coord > y (const double) y field plate coord < meanRA (double *) mean RA < meanDec (double *) mean Declination > equinox (const double) Equinox * Returned value: * See Also: * Support: Scott Smedley, AAO *- * History: 10-Jun-2004 - SS - Initial version. 18-Jun-2004 - SS - Add optional parameter: $equinox. */ void telMdlMean::Pos2Mean ( const double x, const double y, double *meanRA, double *meanDec, const double equinox) const { Pos2App(x, y, meanRA, meanDec); slaAmp(*meanRA, *meanDec, dateMjd, equinox, meanRA, meanDec); } /* *+ t e l M d l M e a n : : P o s 2 M e a n I t e r * Function name: telMdlMean::Pos2MeanIter * Function: Convert field plate coordinates into mean RA/Dec coords. * Description: The telMdlMean::Pos2Mean() member function returns an _approximation_ of the mean RA/Dec coords for a given field plate coordinate. This routine iterates to find a more precise value for the RA/Dec coords. * Call: * Parameters: (">" input, "!" modified, "W" workspace, "<" output) > x (const double) x field plate coord > y (const double) y field plate coord < meanRA (double *) mean RA < meanDec (double *) mean Declination (>) tolerance (double *) accuracy of $ra, $dec coords on field plate. (>) maxIters (double *) maximum number of iterations. > equinox (const double) Equinox * Returned value: * See Also: telMdlMean::Pos2Mean() member function. * Support: Scott Smedley, AAO *- * History: 17-Jun-2004 - SS - Initial version. 18-Jun-2004 - SS - Add optional parameter: $equinox. */ void telMdlMean::Pos2MeanIter ( const double x, const double y, double *meanRA, double *meanDec, const double tolerance, const int maxIters, const double equinox) const { Pos2AppIter(x, y, meanRA, meanDec, tolerance, maxIters); slaAmp(*meanRA, *meanDec, dateMjd, equinox, meanRA, meanDec); } /* *+ t e l M d l M e a n : : E q u * Function name: telMdlMean::DeltaEqu * Function: Determine delta x-y/ra-dec for a given Mean RA/Dec within the field. * Description: This method determines the change in RA & Dec (apparent) indicated by a given change in x &y. The function returns factors such that DeltaRa = DeltaX * drdx + DeltaY * drdy DeltaDec = DeltaY * dddx + DeltaY * dddy This conversion will be local to the specified RA/Dec position (it can't be used in other locations due to field distortion) * Call: * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) meanRa (const double) The RA at which we determine the conversion. Radians. (>) meanDec (const double) The Dec at which we determine the conversion. (>) pmPars (const telMdlPM_PARS &) The proper motion parameters for this position. (<) drdx (double *) Change in RA for a change in X (<) drdy (double *) Change in RA for a change in Y (<) dddx (double *) Change in Dec for a change in X (<) dddy (double *) Change in Dec for a change in Y * Returned value: * See Also: telMdl::DeltaEqu(), telMdlMean::DeltaAltAz(). * Support: Tony Farrell, AAO *- * History: 23-Jan-2008 - THF - Initial version. */ void telMdlMean::DeltaEqu( const double meanRa, const double meanDec, const telMdlPM_PARS &pmPars, double * const drdx, /* Change in RA for a change in x */ double * const drdy, /* Change in RA for a change in y */ double * const dddx, /* Change in Dec for a change in x */ double * const dddy /* Change in Dec for a change in y */ ) const { double ra; double dec; /* * Get apparent RA and Dec values, accounting for proper motion. */ slaMapqk (meanRa, meanDec, pmPars.ra, pmPars.dec, pmPars.parallax, pmPars.radialVel, const_cast(gmap), &ra, &dec); /* * And do the calculation. */ telMdl::DeltaEqu(ra, dec, drdx, drdy, dddx, dddy); } /* *+ t e l M d l M e a n : : E q u * Function name: telMdlMean::DeltaAltAz * Function: Determine delta x-y/az-el for a given Mean RA/Dec within the field. * Description: This method determines the change in Azimuth and Elevation indicated by a given change in x &y. The function returns factors such that DeltaAz = DeltaX * dadx + DeltaY * dady DeltaEl = DeltaY * dedx + DeltaY * dedy This conversion will be local to the specified mean RA/Dec position (it can't be used in other locations due to field distortion). Note - we do mean the specified RA/Dec position, not an Alt/Az position. * Call: * Parameters: (">" input, "!" modified, "W" workspace, "<" output) (>) meanRa (const double) The RA at which we determine the conversion. Radians. (>) meanDec (const double) The Dec at which we determine the conversion. (>) pmPars (const telMdlPM_PARS &) The proper motion parameters for this position. (<) dadx (double *) Change in Azimuth for a change in X (<) dady (double *) Change in Azimuth for a change in Y (<) dedx (double *) Change in RA for a change in X (<) dedy (double *) Change in RA for a change in Y * Returned value: * See Also: telMdl::DeltaAltAz(), telMdlMean::DeltaEqu(). * Support: Tony Farrell, AAO *- * History: 23-Jan-2008 - THF - Initial version. */ void telMdlMean::DeltaAltAz( const double meanRa, const double meanDec, const telMdlPM_PARS &pmPars, double * const dadx, /* Change in Azimuth for a change in x */ double * const dady, /* Change in Azimuth for a change in y */ double * const dedx, /* Change in Elevation for a change in x */ double * const dedy /* Change in Elevation for a change in y */ ) const { double ra; double dec; /* * Get apparent RA and Dec values, accounting for proper motion. */ slaMapqk (meanRa, meanDec, pmPars.ra, pmPars.dec, pmPars.parallax, pmPars.radialVel, const_cast(gmap), &ra, &dec); /* * And do the calculation. */ telMdl::DeltaAltAz(ra, dec, dadx, dady, dedx, dedy); } } // namespace echidna // vim:ts=4 sw=4