Index: /trunk/Ohana/src/libdvo/include/dvo.update.h =================================================================== --- /trunk/Ohana/src/libdvo/include/dvo.update.h (revision 8638) +++ /trunk/Ohana/src/libdvo/include/dvo.update.h (revision 8638) @@ -0,0 +1,399 @@ +# include +# include +# include + +# ifndef DVO_H +# define DVO_H + +/* + # define ELIXIR 1 + # define PANSTARRS 0 + # define LONEOS 0 +*/ + +/*** named data values (convert all to enums?) ***/ + +/* DVO table modes */ +enum {DVO_MODE_UNDEF, DVO_MODE_RAW, DVO_MODE_MEF, DVO_MODE_SPLIT, DVO_MODE_MYSQL} DVOTableMode; + +/* DVO table formats */ +enum {DVO_FORMAT_UNDEF, DVO_FORMAT_INTERNAL, DVO_FORMAT_ELIXIR, DVO_FORMAT_LONEOS, DVO_FORMAT_PANSTARRS, DVO_FORMAT_PMTEST} DVOTableFormat; + +/* image data modes in RegImage */ +enum {T_UNDEF = -1, T_NONE, T_OBJECT, T_DARK, T_BIAS, T_FLAT, T_MASK, T_FRINGE, T_SCATTER, T_MODES, T_FRINGEPTS, T_ANY, N_TYPE}; +enum {M_UNDEF = -1, M_NONE, M_MEF, M_SPLIT, M_SINGLE, M_CUBE, M_SLICE, M_MODES, N_MODE}; + +typedef enum { + PROJ_NONE, // undefined + PROJ_ZEA, // zenithal + PROJ_ZPL, // zenithal + PROJ_ARC, // zenithal + PROJ_STG, // zenithal + PROJ_SIN, // zenithal + PROJ_TAN, // zenithal + PROJ_DIS, // zenithal (TAN + polyterms) + PROJ_LIN, // cartesian + PROJ_PLY, // cartesian + PROJ_WRP, // cartesian + PROJ_AIT, // pseudocyl + PROJ_GLS, // pseudocyl + PROJ_PAR, // pseudocyl +} OhanaProjections; + +typedef enum { + PROJ_MODE_NONE, + PROJ_MODE_CARTESIAN, + PROJ_MODE_ZENITHAL, + PROJ_MODE_PSEUDOCYL, +} OhanaProjectionModes; + +/* RegImage.flag values */ +# define IMREG_DIST 0x01 /* image distributed, only imregister-3.0 */ + +/* catalog values to be loaded */ +# define LOAD_AVES 0x01 +# define LOAD_MEAS 0x02 +# define LOAD_MISS 0x04 +# define LOAD_SECF 0x08 +# define LOAD_MEAS_META 0x10 + +/* invalid mag value */ +# define NO_MAG 0x7fff +# define NO_ERR 0xff + +/* photometry code types */ +# define PHOT_PRI 0x01 +# define PHOT_SEC 0x02 +# define PHOT_DEP 0x03 +# define PHOT_REF 0x04 +# define PHOT_ALT 0x05 /* never stored, only for look-ups */ + +/* Image.code values. these are codes to note bad images */ +# define ID_IMAGE_NEW 0x0000 /* no nrphot attempted */ +# define ID_IMAGE_NOCAL 0x0001 /* used within nrphot to mean "don't apply fit" */ +# define ID_IMAGE_POOR 0x0002 /* relphot says image is bad */ +# define ID_IMAGE_SKIP 0x0004 /* external information image is bad */ +# define ID_IMAGE_FEW 0x0008 /* currently too few measurements for good value */ + +/* Measure.flags values */ +# define ID_MEAS_NOCAL 0x0001 +# define ID_MEAS_POOR 0x0002 +# define ID_MEAS_SKIP 0x0004 +# define ID_MEAS_AREA 0x0008 +# define BLEND_IMAGE 0x0100 +# define BLEND_CATALOG 0x0200 +# define BLEND_IMAGE_NEIGHBOR 0x1000 +# define ID_MEAS_TRAIL 0x2000 +# define ID_MEAS_GHOST 0x4000 + +/* some subtle distinctions between the blend flags: + BLEND_IMAGE: the star on an image is matched with more + than one star in the catalog (image has worse seeing than catalog) + BLEND_CATALOG: the star in the catalog is matched with more + than one star on the image (image has better seeing than catalog) + CALIBRATED: relative photometry has been performed on this measurement + BLEND_IMAGE_NEIGHBOR: the star on an image is matched with more + than one star in the catalog, but not in the same catalog file. +*/ + +/* Average.code values */ +# define ID_STAR_FEW 0x0001 /* used within relphot: skip star */ +# define ID_STAR_POOR 0x0002 /* used within relphot: skip star */ +# define ID_PROPER 0x0400 /* star with large proper motion */ +# define ID_TRANSIENT 0x1000 /* is this mutually exclusive with USNO? */ +# define ID_VARIABLE 0x2000 /* not currently set? */ +# define ID_ASTEROID 0x2000 /* identified with an asteroid */ +# define ID_BAD_OBJECT 0x4000 /* if all measurements are bad, set this bit */ +# define ID_MOVING 0x8000 +# define ID_ROCK 0xa000 /* 0x8000 + 0x2000 */ +# define ID_GHOST 0xc001 /* 0x8000 + 0x4000 + 0x0001 */ +# define ID_TRAIL 0xc002 /* 0x8000 + 0x4000 + 0x0002 */ +# define ID_BLEED 0xc003 /* 0x8000 + 0x4000 + 0x0003 */ +# define ID_COSMIC 0xc004 /* 0x8000 + 0x4000 + 0x0004 */ + +/*** general dvo structures (internal use only / not IO) ***/ + +/* FITS DB structure */ +typedef struct { + FILE *f; + char *filename; + int dbstate; + int lockstate; + double timeout; + Header header; + Matrix matrix; + Header theader; + FTable ftable; + VTable vtable; + int mode; /* what data storage mode is used for disk file? */ + int format; /* what data format is used for disk file? */ + int virtual; /* is table in ftable or vtable? */ + int swapped; /* is table in internal byte-order? */ +} FITS_DB; + +/* the basic HST GSC layout corresponds to a depth of 3 */ +# define SKY_DEPTH_HST 3 + +/* SkyRegion : better implementation than GSCRegion */ +typedef struct { + int Nregions; + char **filename; + SkyRegion *regions; +} SkyTable; + +typedef struct { + int Nregions; + char **filename; + SkyRegion **regions; +} SkyList; + +# if (0) +/* structure for data on a catalog region */ +typedef struct { + char filename[256]; + double DEC[2], RA[2]; +} GSCRegion; +# endif + +typedef struct { + int Ncode; + int Nsecfilt; + int hashcode[0x10000]; + int hashNsec[0x10000]; + PhotCode *code; +} PhotCodeData; + +/* a catalog contains this data */ +typedef struct Catalog { + char *filename; /* catalog file */ + FILE *f; /* file descriptor */ + Header header; + + /* data in the catalog file */ + Average *average; + Measure *measure; + Missing *missing; + SecFilt *secfilt; + int Naverage, Nmeasure, Nmissing, Nsecfilt; /* current number of each component */ + int Nave_disk, Nmeas_disk, Nmiss_disk; /* number of component on disk */ + int Nmeas_off; /* dist seq of first loaded data value */ + /* note the different counting for Nsecfilt */ + + /* pointers to split data files */ + struct Catalog *measure_catalog; /* measure catalog data (split) */ + struct Catalog *missing_catalog; /* missing catalog data (split) */ + struct Catalog *secfilt_catalog; /* secfilt catalog data (split) */ + + /* extra catalog information */ + int lockmode; + int catmode; /* storage mode (raw, mef, split, mysql) */ + int catformat; /* storage format (elixir, panstarrs, etc) */ + int catflags; /* choices to be loaded */ + int sorted; /* is measure table average-sorted? */ + + /* pointers for data manipulation */ + int *found; + int *image; + int *mosaic; + float *X; + float *Y; + +} Catalog; + +/*** prototypes ***/ + +/* in gfits_db.c */ +int gfits_db_init PROTO((FITS_DB *db)); +int gfits_db_create PROTO((FITS_DB *db)); +int gfits_db_lock PROTO((FITS_DB *db, char *filename)); +int gfits_db_load PROTO((FITS_DB *db)); +int gfits_db_load_segment PROTO((FITS_DB *db, int start, int Nrows)); +int gfits_db_save PROTO((FITS_DB *db)); +int gfits_db_update PROTO((FITS_DB *db)); +int gfits_db_close PROTO((FITS_DB *db)); +int gfits_db_free PROTO((FITS_DB *db)); + +/* in coords.c, using libautocode/def/coords.d */ +int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords); +int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords); +int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords); +int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords); +int GetCoords (Coords *coords, Header *header); +int PutCoords (Coords *coords, Header *header); +void RegisterMosaic (Coords *coords); +void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch); + +int FindMosaicForImage (Image *images, int Nimages, int entry); +int FindMosaicForImage_TableSearch (Image *images, int Nimages, int entry); +int FindMosaicForImage_MatchSearch (Image *images, int Nimages, int entry); +int BuildChipMatch (Image *images, int Nimages); +void SetImageCorners (double *X, double *Y, Image *image); + +short int putMi (double value); +double getMi (short int value); +void returnMcal (Image *image, double *c); +void assignMcal (Image *image, double *c, int order); +double applyMcal (Image *image, double x, double y); +double findscatter (double *X, double *Y, double *M, double *dM, int N, double *c, int order); + +PhotCode *GetPhotcodebyName (char *name); +PhotCode *GetPhotcodeEquivbyName (char *name); +PhotCode *GetPhotcodebyCode (int code); +PhotCode *GetPhotcodebyNsec (int Nsec); +PhotCode *GetPhotcodeEquivbyCode (int code); +char *GetPhotcodeNamebyCode (int code); + +float PhotInst (Measure *measure); +float PhotCat (Measure *measure); +float PhotSys (Measure *measure, Average *average, SecFilt *secfilt); +float PhotRel (Measure *measure, Average *average, SecFilt *secfilt); +float PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt); + +# if (0) +short iPhotInst (Measure *measure); +short iPhotAbs (Measure *measure); +short iPhotCat (Measure *measure); +short iPhotSys (Measure *measure, Average *average, SecFilt *secfilt); +short iPhotRel (Measure *measure, Average *average, SecFilt *secfilt); +short iPhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +short iPhotAve (PhotCode *code, Average *average, SecFilt *secfilt); +short iPhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +short iPhotXm (PhotCode *code, Average *average, SecFilt *secfilt); +short iPhotdM (PhotCode *code, Average *average, SecFilt *secfilt); +# endif + +float iPhotColor (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color); + +int LoadPhotcodes (char *filename); +int GetPhotcodeCodebyName (char *name); +int GetPhotcodeEquivCodebyName (char *name); +int GetPhotcodeEquivCodebyCode (int code); +int GetPhotcodeNsec (int code); +int GetPhotcodeNsecfilt (); +void SetZeroPoint (double ZP); +int *GetPhotcodeEquivList (int code, int *nlist); +void ParseColorTerms (char *terms, float *X, int *N); + +int get_image_type (char *name); +char *get_type_name (int type); +int get_image_mode (char *name); +char *get_mode_name (int mode); + +/* these functions refer to the DVO structures defined in the includes above + * they are being replaced with autocode entries (drop when totally autocoded) + */ + +int Fread (void *ptr, int size, int nitems, FILE *f, char *type); +int Fwrite (void *ptr, int size, int nitems, FILE *f, char *type); +int ByteSwap (char *ptr, int size, int nitems, char *type); +int ConvertStruct (char *buffer, int size, int Nbytes, char *type); + +/** dvo_catalog APIs */ +void dvo_catalog_init (Catalog *catalog, int complete); +void dvo_catalog_create (SkyRegion *region, Catalog *catalog); +void dvo_catalog_free (Catalog *catalog); +int dvo_catalog_check (Catalog *catalog, int Nsecfilt, int extend); +int dvo_catalog_lock (Catalog *catalog, int lockmode); +int dvo_catalog_unlock (Catalog *catalog); +int dvo_catalog_load (Catalog *catalog, int VERBOSE); +int dvo_catalog_open (Catalog *catalog, SkyRegion *region, int VERBOSE, char *iomode); +int dvo_catalog_save (Catalog *catalog, char VERBOSE); +int dvo_catalog_update (Catalog *catalog, char VERBOSE); +int dvo_catalog_catformat (char *catformat); +int dvo_catalog_catmode (char *catmode); +void dvo_catalog_test (Catalog *catalog, int halt); + +/* catmode-specific APIs */ +int dvo_catalog_load_raw (Catalog *catalog, int VERBOSE); +int dvo_catalog_save_raw (Catalog *catalog, char VERBOSE); +int dvo_catalog_load_mef (Catalog *catalog, int VERBOSE); +int dvo_catalog_save_mef (Catalog *catalog, char VERBOSE); +int dvo_catalog_load_split (Catalog *catalog, int VERBOSE); +int dvo_catalog_save_split (Catalog *catalog, char VERBOSE); +int dvo_catalog_update_split (Catalog *catalog, char VERBOSE); + +/*** conversion functions / I/O conversions ***/ +Average *ReadRawAverage (FILE *f, int Naverage, int format); +Measure *ReadRawMeasure (FILE *f, int Nmeasure, int format); +SecFilt *ReadRawSecFilt (FILE *f, int Nsecfilt, int format); +int WriteRawAverage (FILE *f, Average *average, int Naverage, int format); +int WriteRawMeasure (FILE *f, Measure *measure, int Nmeasure, int format); +int WriteRawSecFilt (FILE *f, SecFilt *secfilt, int Nsecfilt, int format); + +Average *FtableToAverage (FTable *ftable, int *Naverage, int *format); +Average *AverageLoneosToInternal (AverageLoneos *in, int Nvalues); +Average *AverageElixirToInternal (AverageElixir *in, int Nvalues); +Average *AveragePanstarrsToInternal (AveragePanstarrs *in, int Nvalues); +Average *AveragePMtestToInternal (AveragePMtest *in, int Nvalues); +AverageLoneos *AverageInternalToLoneos (Average *in, int Nvalues); +AverageElixir *AverageInternalToElixir (Average *in, int Nvalues); +AveragePanstarrs *AverageInternalToPanstarrs (Average *in, int Nvalues); +AveragePMtest *AverageInternalToPMtest (Average *in, int Nvalues); + +Measure *FtableToMeasure (FTable *ftable, int *Nmeasure, int *format); +Measure *MeasureLoneosToInternal (MeasureLoneos *in, int Nvalues); +Measure *MeasureElixirToInternal (MeasureElixir *in, int Nvalues); +Measure *MeasurePanstarrsToInternal (MeasurePanstarrs *in, int Nvalues); +MeasureLoneos *MeasureInternalToLoneos (Measure *in, int Nvalues); +MeasureElixir *MeasureInternalToElixir (Measure *in, int Nvalues); +MeasurePanstarrs *MeasureInternalToPanstarrs (Measure *in, int Nvalues); + +SecFilt *FtableToSecFilt (FTable *ftable, int *Nsecfilt, int *format); +SecFilt *SecFiltLoneosToInternal (SecFiltLoneos *in, int Nvalues); +SecFilt *SecFiltElixirToInternal (SecFiltElixir *in, int Nvalues); +SecFilt *SecFiltPanstarrsToInternal (SecFiltPanstarrs *in, int Nvalues); +SecFiltLoneos *SecFiltInternalToLoneos (SecFilt *in, int Nvalues); +SecFiltElixir *SecFiltInternalToElixir (SecFilt *in, int Nvalues); +SecFiltPanstarrs *SecFiltInternalToPanstarrs (SecFilt *in, int Nvalues); + +int AverageToFtable (FTable *ftable, Average *average, int Naverage, int format); +int MeasureToFtable (FTable *ftable, Measure *measure, int Nmeasure, int format); +int SecFiltToFtable (FTable *ftable, SecFilt *secfilt, int Nsecfilt, int format); + +/*** DVO image db I/O Functions ***/ +int dvo_image_lock (FITS_DB *db, char *filename, double timeout, int lockstate); +int dvo_image_unlock (FITS_DB *db); +int dvo_image_load (FITS_DB *db, int VERBOSE, int FORCE_READ); +int dvo_image_save (FITS_DB *db, int VERBOSE); +int dvo_image_update (FITS_DB *db, int VERBOSE); +int dvo_image_load_raw (FITS_DB *db, int VERBOSE, int FORCE_READ); +int dvo_image_update_raw (FITS_DB *db, int VERBOSE); +int dvo_image_save_raw (FITS_DB *db, int VERBOSE); +int dvo_image_addrows (FITS_DB *db, Image *new, int Nnew); +void dvo_image_create (FITS_DB *db, double ZeroPoint); + +int FtableToImage (FTable *ftable, Header *theader, int *format); +int ImageToFtable (FTable *ftable, Header *theader, int format); +int ImageToVtable (VTable *vtable, Header *theader, int format); +Image *ImageElixirToInternal (ImageElixir *in, int Nvalues); +ImageElixir *ImageInternalToElixir (Image *in, int Nvalues); +Image *ImageLoneosToInternal (ImageLoneos *in, int Nvalues); +ImageLoneos *ImageInternalToLoneos (Image *in, int Nvalues); +Image *ImagePanstarrsToInternal (ImagePanstarrs *in, int Nvalues); +ImagePanstarrs *ImageInternalToPanstarrs (Image *in, int Nvalues); + +/* skyregion APIs */ +int SkyTableSave PROTO((SkyTable *table, char *filename)); +SkyTable *SkyTableLoad PROTO((char *filename, int VERBOSE)); +SkyTable *SkyTableFromGSC PROTO((char *filename, int depth, int VERBOSE)); +SkyTable *SkyTableLoadOptimal PROTO((char *catdir, char *SKYFILE, char *GSCFILE, int depth, int VERBOSE)); +int SkyTableSetDepth PROTO((SkyTable *sky, int depth)); +SkyList *SkyRegionByPoint PROTO((SkyTable *table, int depth, double ra, double dec)); +SkyList *SkyListByPoint PROTO((SkyTable *table, double ra, double dec)); +SkyList *SkyListByRadius PROTO((SkyTable *table, int depth, double RA, double DEC, double radius)); +SkyList *SkyListByPatch PROTO((SkyTable *table, int depth, SkyRegion *patch)); +SkyList *SkyListByName PROTO((SkyTable *table, char *name)); +SkyList *SkyListByImage PROTO((SkyTable *table, int depth, Image *image)); +SkyList *SkyListByBounds PROTO((SkyTable *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax)); +SkyList *SkyListChildrenByBounds PROTO((SkyTable *table, int No, int depth, double Rmin, double Rmax, double Dmin, double Dmax)); +int SkyListFree PROTO((SkyList *list, int ELEMENTS)); +int SkyTableFree PROTO((SkyTable *table)); +int SkyListSetFilenames PROTO((SkyList *list, char *path, char *ext)); +int SkyTableSetFilenames PROTO((SkyTable *sky, char *path, char *ext)); + +# endif Index: /trunk/Ohana/src/libdvo/src/coordops.update.c =================================================================== --- /trunk/Ohana/src/libdvo/src/coordops.update.c (revision 8638) +++ /trunk/Ohana/src/libdvo/src/coordops.update.c (revision 8638) @@ -0,0 +1,657 @@ +# include + +static Coords *mosaic = NULL; + +void RegisterMosaic (Coords *coords) { + mosaic = coords; +} + +int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) { + + int Zenith1, Zenith2, Zenithal, Polynomial, Cartesian, PseudoCyl; + char *type; + double L, M, X, Y, T, Z, Z2; + double R, sphi, cphi, stht, ctht; + double alpha, delta, salp, calp, sdel, sdp, cdp; + + *ra = 0; + *dec = 0; + stht = ctht = 1; + + proj = GetProjection (coords[0].ctype); + mode = GetProjectionMode (proj); + if (proj == PROJ_NONE) return (FALSE); + if (proj == PROJ_MODE_NONE) return (FALSE); + + /** convert pixel coordinates to cartesian system **/ + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + + /** extra polynomial terms **/ + if (coords[0].Npolyterms > 1) { + L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0]; + M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1]; + } + if (coords[0].Npolyterms > 2) { + L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0]; + M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1]; + } + + /** Locally Cartesian Projections **/ + if (mode == PROJ_MODE_CARTESIAN) { + *ra = L + coords[0].crval1; + *dec = M + coords[0].crval2; + + /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */ + if (proj == PROJ_WRP) { + if (mosaic == NULL) return (FALSE); + XY_to_RD (ra, dec, L + coords[0].crval1, M + coords[0].crval2, mosaic); + } + return (TRUE); + } + + /** Zenithal Projections **/ + if (mode == PROJ_MODE_ZENITHAL) { + R = hypot (L,M); + if ((L == 0) && (M == 0)) { + sphi = 0; + cphi = 1; + } else { + sphi = L / R; + cphi = -M / R; + } + + switch (proj) { + case PROJ_TAN: + if (R == 0) { + stht = 1.0; + ctht = 0.0; + } else { + T = DEG_RAD / R; + stht = T / sqrt ( 1.0 + T*T); + ctht = 1.0 / sqrt ( 1.0 + T*T); + } + break; + case PROJ_STG: + stht = (4 - RAD_DEG*R) / (4 + RAD_DEC*R); + ctht = sqrt (1 - stht*stht); + break; + case PROJ_SIN: + ctht = RAD_DEG * R; + stht = sqrt (1 - ctht*ctht); + break; + case PROJ_ZEA: + case PROJ_ZPL: + stht = 1 - 0.5*SQ(R*RAD_DEG); + ctht = sqrt (1 - stht*stht); + break; + } + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + + sdel = stht*sdp - ctht*cphi*cdp; + salp = ctht*sphi; + calp = stht*cdp + ctht*cphi*sdp; + alpha = atan2 (salp, calp); + delta = asin (sdel); + + *ra = DEG_RAD*alpha + coords[0].crval1; + *dec = DEG_RAD*delta; + + /* rationalize ra range 0 - 360.0 */ + while (*ra < 0.0) *ra += 360.0; + while (*ra > 360.0) *ra -= 360.0; + + return (TRUE); + } + + /**** Other Conventional Projections ****/ + if (mode == PROJ_MODE_ZENITHAL) { + switch (proj) { + case PROJ_AIT: + Z2 = (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M)); + if (Z2 < 0) return (FALSE); + Z = sqrt (Z2); + alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*Z2 - 1.0); + delta = DEG_RAD * asin (RAD_DEG*M*Z); + break; + + case PROJ_GLS: + /* L,M in degrees, alpha,delta in degrees */ + alpha = L / cos (RAD_DEG * M); + delta = M; + break; + case PROJ_PAR: + /* L,M in degrees, alpha,delta in degrees */ + alpha = L / (1.0 - SQ(2.0*M/180)); + delta = 3 * DEG_RAD * asin (M/180.0); + break; + } + *ra = alpha + coords[0].crval1; + *dec = delta + coords[0].crval2; + + /* rationalize ra range 0 - 360.0 */ + while (*ra < 0.0) *ra += 360.0; + while (*ra > 360.0) *ra -= 360.0; + + return (TRUE); + } + return (FALSE); +} + +int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords) { + + char *type; + int i, status, Polynomial, Zenith1, Zenith2, Zenithal, Cartesian, PseudoCyl; + double phi, theta; + double determ; + double X, Y, L, M, Lo, Mo, dL, dM; + double sphi, cphi, stht; + double salp, calp, sdel, cdel, sdp, cdp; + double P, A, Rc; + + status = TRUE; + *x = 0; + *y = 0; + type = &coords[0].ctype[4]; + L = M = 0; + + proj = GetProjection (coords[0].ctype); + mode = GetProjectionMode (proj); + if (proj == PROJ_NONE) return (FALSE); + if (proj == PROJ_MODE_NONE) return (FALSE); + + /**** Locally Cartesian Projections ****/ + if (mode == PROJ_MODE_CARTESIAN) { + if (proj == PROJ_WRP) { + if (mosaic == NULL) return (FALSE); + RD_to_XY (&Lo, &Mo, ra, dec, mosaic); + L = (Lo - coords[0].crval1); + M = (Mo - coords[0].crval2); + } else { + L = (ra - coords[0].crval1); + M = (dec - coords[0].crval2); + } + } + + /**** Zenithal Projections ****/ + if (mode == PROJ_MODE_ZENITHAL) { + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + salp = sin(RAD_DEG*(ra - coords[0].crval1)); + calp = cos(RAD_DEG*(ra - coords[0].crval1)); + sdel = sin(RAD_DEG*dec); + cdel = cos(RAD_DEG*dec); + + stht = sdel*sdp + cdel*cdp*calp; /* sin(theta) */ + sphi = cdel*salp; /* = cos(theta)*sin(phi) */ + cphi = cdel*sdp*calp - sdel*cdp; /* = cos(theta)*cos(phi) */ + if (stht < 0) status = FALSE; + + switch (proj) { + case PROJ_TAN: + case PROJ_DIS: + L = DEG_RAD * sphi / stht; + M = -DEG_RAD * cphi / stht; + break; + case PROJ_SIN: + L = DEG_RAD * sphi; + M = -DEG_RAD * cphi; + break; + case PROJ_ZEA: + case PROJ_ZPL: + Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht); + L = Rc * sphi; + M = -Rc * cphi; + status = TRUE; + break; + } + } + + /**** Other Standard Projections ****/ + if (PseudoCyl) { + if (!strcmp(type, "-AIT")) { + phi = RAD_DEG*(ra - coords[0].crval1); + theta = RAD_DEG*(dec - coords[0].crval2); + P = 1.0 + cos (theta) * cos (0.5*phi); + if (P != 0.0) { + A = DEG_RAD * sqrt (2.0 / P); + L = 2.0 * A * cos (theta) * sin (0.5*phi); + M = A * sin (theta); + } else { + L = 0.0; + M = 0.0; + } + } + if (!strcmp(type, "-GLS")) { + phi = ra - coords[0].crval1; + theta = dec - coords[0].crval2; + L = phi * cos(RAD_DEG * theta); + M = theta; + } + if (!strcmp(type, "-PAR")) { + phi = ra - coords[0].crval1; + theta = dec - coords[0].crval2; + L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1); + M = 180.0 * sin (RAD_DEG*theta/3.0); + } + } + + /* convert L,M to X,Y */ + determ = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1); + X = determ * (coords[0].pc2_2*L - coords[0].pc1_2*M); + Y = determ * (coords[0].pc1_1*M - coords[0].pc2_1*L); + + /** extra polynomial terms **/ + if (coords[0].Npolyterms > 1) { + for (i = 0; i < 10; i++) { + Lo = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + Mo = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + if (coords[0].Npolyterms > 1) { + Lo += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0]; + Mo += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1]; + } + if (coords[0].Npolyterms > 2) { + Lo += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0]; + Mo += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1]; + } + dL = (L - Lo); + dM = (M - Mo); + // fprintf (stderr, "%d: %f,%f : %f,%f : %f,%f : %f,%f\n", i, L, M, X, Y, Lo, Mo, dL, dM); + + X += determ * (coords[0].pc2_2*dL - coords[0].pc1_2*dM); + Y += determ * (coords[0].pc1_1*dM - coords[0].pc2_1*dM); + } + } + /* check for correct size (iterate?) */ + + *x = X / coords[0].cdelt1 + coords[0].crpix1; + *y = Y / coords[0].cdelt2 + coords[0].crpix2; + + return (status); + +} + +int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords) { + + int status; + double tmpx, tmpy; + + status = RD_to_XY (&tmpx, &tmpy, ra, dec, coords); + *x = tmpx; + *y = tmpy; + + return (status); + +} + +int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords) { + + int status; + double tmpr, tmpd; + + status = XY_to_RD (&tmpr, &tmpd, x, y, coords); + *ra = tmpr; + *dec = tmpd; + + return (status); + +} + +int GetCoords (Coords *coords, Header *header) { + + int status, itmp, Polynomial, Polyterm; + double Lambda, rotate, scale; + double equinox; + char *ctype; + + rotate = 0.0; + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + coords[0].Npolyterms = 1; + strcpy (coords[0].ctype, "NONE"); + + status = FALSE; + if (gfits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) { + status = gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); + status &= gfits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1); + status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); + status &= gfits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2); + + if (gfits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1)) { + status &= gfits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2); + if (gfits_scan (header, "CROTA2", "%lf", 1, &rotate)) { + Lambda = coords[0].cdelt2 / coords[0].cdelt1; + coords[0].pc1_1 = cos(rotate*RAD_DEG); + coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda; + coords[0].pc2_1 = sin(rotate*RAD_DEG) / Lambda; + coords[0].pc2_2 = cos(rotate*RAD_DEG); + } + if (gfits_scan (header, "PC001001", "%f", 1, &coords[0].pc1_1)) { + status &= gfits_scan (header, "PC001002", "%f", 1, &coords[0].pc1_2); + status &= gfits_scan (header, "PC002001", "%f", 1, &coords[0].pc2_1); + status &= gfits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2); + } + + /* set NPLYTERM based on header. if NPLYTERM is missing, it should have a + value of 0, unless the projection type is one of PLY, DIS, WRP, in which + case it should be set to 3 */ + ctype = &coords[0].ctype[4]; + Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP"); + Polyterm = gfits_scan (header, "NPLYTERM", "%d", 1, &itmp); + + coords[0].Npolyterms = 0; + if (Polynomial && !Polyterm) coords[0].Npolyterms = 3; + if (Polyterm) coords[0].Npolyterms = itmp; + + switch (coords[0].Npolyterms) { + case 3: + status &= gfits_scan (header, "PCA1X3Y0", "%f", 1, &coords[0].polyterms[3][0]); + status &= gfits_scan (header, "PCA1X2Y1", "%f", 1, &coords[0].polyterms[4][0]); + status &= gfits_scan (header, "PCA1X1Y2", "%f", 1, &coords[0].polyterms[5][0]); + status &= gfits_scan (header, "PCA1X0Y3", "%f", 1, &coords[0].polyterms[6][0]); + status &= gfits_scan (header, "PCA2X3Y0", "%f", 1, &coords[0].polyterms[3][1]); + status &= gfits_scan (header, "PCA2X2Y1", "%f", 1, &coords[0].polyterms[4][1]); + status &= gfits_scan (header, "PCA2X1Y2", "%f", 1, &coords[0].polyterms[5][1]); + status &= gfits_scan (header, "PCA2X0Y3", "%f", 1, &coords[0].polyterms[6][1]); + case 2: + status &= gfits_scan (header, "PCA1X2Y0", "%f", 1, &coords[0].polyterms[0][0]); + status &= gfits_scan (header, "PCA1X1Y1", "%f", 1, &coords[0].polyterms[1][0]); + status &= gfits_scan (header, "PCA1X0Y2", "%f", 1, &coords[0].polyterms[2][0]); + status &= gfits_scan (header, "PCA2X2Y0", "%f", 1, &coords[0].polyterms[0][1]); + status &= gfits_scan (header, "PCA2X1Y1", "%f", 1, &coords[0].polyterms[1][1]); + status &= gfits_scan (header, "PCA2X0Y2", "%f", 1, &coords[0].polyterms[2][1]); + case 0: + case 1: + break; + } + } else { + if (gfits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1)) { + status &= gfits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2); + status &= gfits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1); + status &= gfits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2); + /* renormalize */ + scale = hypot (coords[0].pc1_1, coords[0].pc1_2); + coords[0].cdelt1 = coords[0].cdelt2 = scale; + coords[0].pc1_1 /= scale; + coords[0].pc1_2 /= scale; + coords[0].pc2_1 /= scale; + coords[0].pc2_2 /= scale; + } else { + status = FALSE; + } + } + } else { + /* some of my thesis data uses this simple linear model - convert on read? */ + if (gfits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) { + status = gfits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1); + status &= gfits_scan (header, "RA_Y", "%f", 1, &coords[0].pc1_2); + status &= gfits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2); + status &= gfits_scan (header, "DEC_X", "%f", 1, &coords[0].pc2_1); + status &= gfits_scan (header, "DEC_Y", "%f", 1, &coords[0].pc2_2); + coords[0].crpix1 = coords[0].crpix2 = 0.0; + coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + strcpy (coords[0].ctype, "GENE"); + } + } + if (status) { + if (!gfits_scan (header, "EQUINOX", "%lf", 1, &equinox)) { + if (!gfits_scan (header, "EPOCH", "%lf", 1, &equinox)) { + equinox = 2000.0; + } + } + if (fabs (equinox - 2000.0) > 0.1) { + coords_precess (&coords[0].crval1, &coords[0].crval2, equinox, 2000.0); + } + } + if (!status) { + fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype); + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + strcpy (coords[0].ctype, "NONE"); + } + return (status); +} + +int PutCoords (Coords *coords, Header *header) { + + int OldAIPS; + char csys[16], ctype[16]; + double rotate, Lambda; + + /* modifications to the ctype? */ + OldAIPS = FALSE; + gfits_modify (header, "CTYPE2", "%s", 1, coords[0].ctype); + if (!strcmp(coords[0].ctype, "MM")) { + gfits_modify (header, "CTYPE1", "%s", 1, "LL"); + OldAIPS = TRUE; + } else { + strcpy (csys, "NONE"); + if (!strncmp (coords[0].ctype, "DEC-", 4)) strcpy (csys, "RA--"); + if (!strncmp (coords[0].ctype, "GLAT", 4)) strcpy (csys, "GLON"); + if (!strncmp (coords[0].ctype, "ELAT", 4)) strcpy (csys, "ELON"); + if (!strncmp (coords[0].ctype, "HLAT", 4)) strcpy (csys, "HLON"); + if (!strncmp (coords[0].ctype, "SLAT", 4)) strcpy (csys, "SLON"); + if (!strcmp (csys, "NONE")) return (FALSE); + sprintf (ctype, "%s-%s", csys, &coords[0].ctype[5]); + gfits_modify (header, "CTYPE1", "%s", 1, ctype); + } + + gfits_modify (header, "CDELT1", "%le", 1, coords[0].cdelt1); + gfits_modify (header, "CDELT2", "%le", 1, coords[0].cdelt2); + gfits_modify (header, "CRVAL1", "%lf", 1, coords[0].crval1); + gfits_modify (header, "CRVAL2", "%lf", 1, coords[0].crval2); + gfits_modify (header, "CRPIX1", "%lf", 1, coords[0].crpix1); + gfits_modify (header, "CRPIX2", "%lf", 1, coords[0].crpix2); + + if (OldAIPS) { + Lambda = coords[0].cdelt2 / coords[0].cdelt1; + rotate = DEG_RAD*atan2 (coords[0].pc2_1*Lambda, coords[0].pc1_1); + gfits_modify (header, "CROTA1", "%f", 1, rotate); + gfits_modify (header, "CROTA2", "%f", 1, rotate); + return (TRUE); + } + + gfits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1); + gfits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2); + gfits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1); + gfits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2); + gfits_modify (header, "NPLYTERM", "%d", 1, coords[0].Npolyterms); + + /* RA Terms */ + if (coords[0].Npolyterms > 1) { + gfits_modify (header, "PCA1X2Y0", "%le", 1, coords[0].polyterms[0][0]); /* polyterms[0]); */ + gfits_modify (header, "PCA1X1Y1", "%le", 1, coords[0].polyterms[1][0]); /* polyterms[1]); */ + gfits_modify (header, "PCA1X0Y2", "%le", 1, coords[0].polyterms[2][0]); /* polyterms[2]); */ + } + if (coords[0].Npolyterms > 2) { + gfits_modify (header, "PCA1X3Y0", "%le", 1, coords[0].polyterms[3][0]); /* polyterms[3]); */ + gfits_modify (header, "PCA1X2Y1", "%le", 1, coords[0].polyterms[4][0]); /* polyterms[4]); */ + gfits_modify (header, "PCA1X1Y2", "%le", 1, coords[0].polyterms[5][0]); /* polyterms[5]); */ + gfits_modify (header, "PCA1X0Y3", "%le", 1, coords[0].polyterms[6][0]); /* polyterms[6]); */ + } + + /* Dec Terms */ + if (coords[0].Npolyterms > 1) { + gfits_modify (header, "PCA2X2Y0", "%le", 1, coords[0].polyterms[0][1]); /* polyterms[7]); */ + gfits_modify (header, "PCA2X1Y1", "%le", 1, coords[0].polyterms[1][1]); /* polyterms[8]); */ + gfits_modify (header, "PCA2X0Y2", "%le", 1, coords[0].polyterms[2][1]); /* polyterms[9]); */ + } + if (coords[0].Npolyterms > 2) { + gfits_modify (header, "PCA2X3Y0", "%le", 1, coords[0].polyterms[3][1]); /* polyterms[10]); */ + gfits_modify (header, "PCA2X2Y1", "%le", 1, coords[0].polyterms[4][1]); /* polyterms[11]); */ + gfits_modify (header, "PCA2X1Y2", "%le", 1, coords[0].polyterms[5][1]); /* polyterms[12]); */ + gfits_modify (header, "PCA2X0Y3", "%le", 1, coords[0].polyterms[6][1]); /* polyterms[13]); */ + } + return (TRUE); +} + +void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch) { + + double T; + double A, D, RA, DEC, zeta, z, theta; + double SA, CA, SD, CD; + + T = (out_epoch - in_epoch) / 100.0; + + zeta = RAD_DEG*(0.6406161*T + 0.0000839*T*T + 0.0000050*T*T*T); + theta = RAD_DEG*(0.5567530*T - 0.0001185*T*T - 0.0000116*T*T*T); + z = 0.6406161*T + 0.0003041*T*T + 0.0000051*T*T*T; + + A = *ra; + D = *dec; + SD = cos(RAD_DEG*A + zeta)*sin(theta)*cos(RAD_DEG*D) + cos(theta)*sin(RAD_DEG*D); + CD = sqrt (1 - SD*SD); + SA = sin(RAD_DEG*A + zeta)*cos(RAD_DEG*D)/CD; + CA = (cos(RAD_DEG*A + zeta)*cos(theta)*cos(RAD_DEG*D) - sin(theta)*sin(RAD_DEG*D))/CD; + + DEC = DEG_RAD*asin(SD); + RA = DEG_RAD*atan2(SA, CA) + z; + + if (RA < 0) + RA += 360; + + *ra = RA; + *dec = DEC; +} + +/* -PLY projection is an extrapolation of the -TAN projection. + In addition to the usual linear terms of CRPIXi, PC00i00j, there are + higher order polynomial terms (up to 3rd order): + Axis 1 terms: + PCA1X2Y0 = coords.polyterm[0][0] = x^2 + PCA1X1Y1 = coords.polyterm[1][0] = xy + PCA1X0Y2 = coords.polyterm[2][0] = y^2 + PCA1X3Y0 = coords.polyterm[3][0] = x^3 + PCA1X2Y1 = coords.polyterm[4][0] = x^2 y + PCA1X1Y2 = coords.polyterm[5][0] = x y^2 + PCA1X0Y3 = coords.polyterm[6][0] = y^2 + Axis 2 terms: + PCA2X2Y0 = coords.polyterm[0][1] = x^2 + PCA2X1Y1 = coords.polyterm[1][1] = xy + PCA2X0Y2 = coords.polyterm[2][1] = y^2 + PCA2X3Y0 = coords.polyterm[3][1] = x^3 + PCA2X2Y1 = coords.polyterm[4][1] = x^2 y + PCA2X1Y2 = coords.polyterm[5][1] = x y^2 + PCA2X0Y3 = coords.polyterm[6][1] = y^2 +*/ + +# if (0) + + /** convert pixel coordinates to cartesian system **/ + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + if (Polynomi) { + if (coords[0].Npolyterms > 2) { + X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]); + Y += coords[0].cdelt2*(x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1]); + } + if (coords[0].Npolyterms > 2) { + X += coords[0].cdelt1*(x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]); + Y += coords[0].cdelt2*(x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]); + } + } + + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); +/** this code is the old method used for higher order terms. they + are essentially 6th order, with weird coupled terms. + I don't think any real data used these terms, but they should + be re-calculated, I would think +**/ + +# endif + +/* + +Projections and Transformations + +The Coords structure is used to represent the standard FITS header representations of the WCS terms. +The FITS WCS encapsulates several concepts in one system: coordinate frame, projection, and +transformations. The ctype variable defines both the frame (ie, RA/DEC, GLON/GLAT, etc) and the +projection (ie, AIT, SIN, TAN, etc). The associated terms (PC00i00j) define transformations from +the projection system to another linear coordinate system. I have extended the basic WCS +transformation terms to include higher-order polynomial terms. The presence of higher-order terms +is indicated in the header by the NPLYTERM keyword, with a value greater than 1 (a value of 0 or 1 +implies no higher-order coefficients). The coefficients have keywords of the form PCAnXiYj where +'n' is the axis corresponding to the dependent variable, 'i' is the order of the X component and 'j' +is the order of the Y component. A value of 2 indicates that the second-order coefficients are +defined (PCAnX2Y0, PCAnX1Y1, PCAnX0Y2); A value of 3 indicates that the third-order coefficients are +also defined (PCAnX3Y0, PCAnX2Y1, PCAnX1Y2, PCAnX0Y3). Some headers in the past were generated +without the NPLYTERM keyword. a value of PLY, DIS, or WRP for the projection without a +corresponding value for NPLYTERM implies that the value should be interpreted as '3'. + +*/ + + /* PLY is equiv to LIN with higher order terms + ZPL is equiv to ZEA with higher order terms + DIS is equiv to TAN with higher order terms + WRP is equiv to PLY, with implied mosaic */ + +int GetProjection (char *ctype) { + if (!strcmp(&ctype[4], "-ZEA")) return PROJ_ZEA; + if (!strcmp(&ctype[4], "-ZPL")) return PROJ_ZPL; + if (!strcmp(&ctype[4], "-ARC")) return PROJ_ARC; + if (!strcmp(&ctype[4], "-STG")) return PROJ_STG; + if (!strcmp(&ctype[4], "-SIN")) return PROJ_SIN; + if (!strcmp(&ctype[0], "MM")) return PROJ_SIN; // note ctype[0] + if (!strcmp(&ctype[4], "-TAN")) return PROJ_TAN; + if (!strcmp(&ctype[4], "-DIS")) return PROJ_DIS; + if (!strcmp(&ctype[4], "-LIN")) return PROJ_LIN; + if (!strcmp(&ctype[0], "GENE")) return PROJ_LIN; // note ctype[0] + if (!strcmp(&ctype[4], "-PLY")) return PROJ_PLY; + if (!strcmp(&ctype[4], "-WRP")) return PROJ_WRP; + if (!strcmp(&ctype[4], "-AIT")) return PROJ_AIT; + if (!strcmp(&ctype[4], "-GLS")) return PROJ_GLS; + if (!strcmp(&ctype[4], "-PAR")) return PROJ_PAR; + return PROJ_NONE; +} + +int SetProjection (char *ctype, int proj) { + switch (proj) { + case PROJ_ZEA: strcpy(&ctype[4], "-ZEA") return TRUE; + case PROJ_ZPL: strcpy(&ctype[4], "-ZPL") return TRUE; + case PROJ_ARC: strcpy(&ctype[4], "-ARC") return TRUE; + case PROJ_STG: strcpy(&ctype[4], "-STG") return TRUE; + case PROJ_SIN: strcpy(&ctype[4], "-SIN") return TRUE; + case PROJ_TAN: strcpy(&ctype[4], "-TAN") return TRUE; + case PROJ_DIS: strcpy(&ctype[4], "-DIS") return TRUE; + case PROJ_LIN: strcpy(&ctype[4], "-LIN") return TRUE; + case PROJ_PLY: strcpy(&ctype[4], "-PLY") return TRUE; + case PROJ_WRP: strcpy(&ctype[4], "-WRP") return TRUE; + case PROJ_AIT: strcpy(&ctype[4], "-AIT") return TRUE; + case PROJ_GLS: strcpy(&ctype[4], "-GLS") return TRUE; + case PROJ_PAR: strcpy(&ctype[4], "-PAR") return TRUE; + } + return FALSE; +} + +int GetProjectionMode (int proj) { + switch (proj) { + case PROJ_ZEA: + case PROJ_ZPL: + case PROJ_ARC: + case PROJ_STG: + case PROJ_SIN: + case PROJ_TAN: + case PROJ_DIS: + return PROJ_MODE_ZENITHAL; + case PROJ_LIN: + case PROJ_PLY: + case PROJ_WRP: + return PROJ_MODE_CARTESIAN; + case PROJ_AIT: + case PROJ_GLS: + case PROJ_PAR: + return PROJ_MODE_PSEUDOCLY; + } + return PROJ_MODE_NONE; +} +