IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10823


Ignore:
Timestamp:
Dec 22, 2006, 11:19:47 AM (20 years ago)
Author:
magnier
Message:

added a bunch of needed projection varients

Location:
trunk/psLib/src/astro
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/psCoord.c

    r10609 r10823  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.127 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2006-12-10 17:29:29 $
     12*  @version $Revision: 1.128 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2006-12-22 21:19:47 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2626#include "psLogMsg.h"
    2727#include "psTrace.h"
     28#include "psAbort.h"
    2829
    2930#include "psMatrix.h"
    … …  
    347348}
    348349
    349 
     350typedef enum {
     351    PS_PROJ_NONE,
     352    PS_PROJ_CARTESIAN,
     353    PS_PROJ_ZENITHAL,
     354    PS_PROJ_PSEUDOCYL,
     355}psProjectionClass;
    350356
    351357psPlane* psProject(
    … …  
    364370    psF64 sinDp, cosDp, sinAlpha, cosAlpha, sinDelta, cosDelta, zeta;
    365371
    366     bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN);
     372    psProjectionClass class = PS_PROJ_NONE;
     373    switch (projection->type) {
     374    case PS_PROJ_LIN:
     375    case PS_PROJ_PLY:
     376    case PS_PROJ_WRP:
     377        class = PS_PROJ_CARTESIAN;
     378        break;
     379    case PS_PROJ_TAN:
     380    case PS_PROJ_DIS:
     381    case PS_PROJ_SIN:
     382    case PS_PROJ_STG:
     383    case PS_PROJ_ZEA:
     384    case PS_PROJ_ZPL:
     385        class = PS_PROJ_ZENITHAL;
     386        break;
     387    case PS_PROJ_AIT:
     388    case PS_PROJ_PAR:
     389    case PS_PROJ_GLS:
     390    case PS_PROJ_MER:
     391    case PS_PROJ_CAR:
     392        class = PS_PROJ_PSEUDOCYL;
     393        break;
     394
     395    default:
     396        psAbort ("psCoord", "invalid projection");
     397        break;
     398    }
    367399
    368400    // Allocate return value
    … …  
    374406    }
    375407
    376     if (zenithal) {
     408    switch (class) {
     409    case PS_PROJ_CARTESIAN:
     410        out->x = (coord->r - projection->R);
     411        out->y = (coord->d - projection->D);
     412        break;
     413
     414    case PS_PROJ_ZENITHAL:
    377415        sinDp = sin(projection->D);
    378416        cosDp = cos(projection->D);
    … …  
    389427    # else
    390428
    391         sinTheta =  sinDelta*sinDp + cosDelta*cosDp*cosAlpha;
    392     cosPhiCT =  sinDelta*cosDp - cosDelta*sinDp*cosAlpha;
     429        sinTheta =  cosDelta*cosAlpha*cosDp + sinDelta*sinDp;
    393430    sinPhiCT = -cosDelta*sinAlpha;
     431    cosPhiCT = -cosDelta*cosAlpha*sinDp + sinDelta*cosDp;
    394432    # endif
    395 
    396 } else {
    397     phi = coord->r - projection->R;
    398     theta = coord->d - projection->D;
    399 }
    400 
    401 // Perform the specified projection
    402 switch (projection->type) {
    403     case PS_PROJ_TAN:
    404         // Gnomonic projection
    405         out->x = +sinPhiCT / sinTheta;
    406         out->y = -cosPhiCT / sinTheta;
    407         break;
    408     case PS_PROJ_SIN:
    409         // Othrographic projection
    410         out->x = +sinPhiCT;
    411         out->y = -cosPhiCT;
    412         break;
    413     case PS_PROJ_AIT:
    414         // Hammer-Aitoff projection
    415         zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0)));
    416         out->x = 2.0*zeta*cos(theta)*sin(phi/2.0);
    417         out->y = zeta*sin(theta);
    418         break;
    419     case PS_PROJ_PAR:
    420         // Parabolic projection
    421         out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0);
    422         out->y = M_PI*sin(theta/3.0);
     433    // Perform the specified projection
     434    switch (projection->type)
     435        {
     436        case PS_PROJ_TAN:
     437        case PS_PROJ_DIS:
     438            // Gnomonic projection
     439            out->x = +sinPhiCT / sinTheta;
     440            out->y = -cosPhiCT / sinTheta;
     441            break;
     442        case PS_PROJ_SIN:
     443            // Othrographic projection
     444            out->x = +sinPhiCT;
     445            out->y = -cosPhiCT;
     446            break;
     447        case PS_PROJ_STG:
     448            // Othrographic projection
     449            out->x = +2*sinPhiCT / (1 + sinTheta);
     450            out->y = -2*cosPhiCT / (1 + sinTheta);
     451            break;
     452        case PS_PROJ_ZEA:
     453        case PS_PROJ_ZPL:
     454            // Zenithal Equal-Area
     455            zeta = M_SQRT2 / sqrt (1 + sinTheta);
     456            out->x = +zeta * sinPhiCT;
     457            out->y = -zeta * cosPhiCT;
     458            break;
     459        default:
     460            psAbort ("psCoord", "invalid projection");
     461            break;
     462        }
     463        break;
     464
     465    case PS_PROJ_PSEUDOCYL:
     466        phi = coord->r - projection->R;
     467        theta = coord->d - projection->D;
     468        switch (projection->type)
     469        {
     470        case PS_PROJ_AIT:
     471            // Hammer-Aitoff projection
     472            zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0)));
     473            out->x = 2.0*zeta*cos(theta)*sin(phi/2.0);
     474            out->y = zeta*sin(theta);
     475            break;
     476        case PS_PROJ_GLS:
     477            // projection name?
     478            out->x = phi * cos(theta);
     479            out->x = theta;
     480            break;
     481        case PS_PROJ_PAR:
     482            // Parabolic projection
     483            out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0);
     484            out->y = M_PI*sin(theta/3.0);
     485            break;
     486        case PS_PROJ_CAR:
     487        case PS_PROJ_MER:
     488            psAbort ("psCoord", "projection not yet implemented");
     489            break;
     490        default:
     491            psAbort ("psCoord", "invalid projection");
     492            break;
     493        }
     494        break;
     495
    423496    default:
    424         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    425                 _("The projection type, %d, is unknown."),
    426                 projection->type);
    427         psFree(out);
    428         return NULL;
     497        psAbort ("psCoord", "invalid projection");
     498        break;
    429499    }
    430500
    … …  
    446516
    447517    psF64 rho      = 0.0;
     518    psF64 rho2      = 0.0;
    448519    psF64 sinTheta = 0.0;
    449520    psF64 cosTheta = 0.0;
    … …  
    453524    psF64  theta = 0.0;
    454525    psF64  phi   = 0.0;
     526    psF64  x, y, R;
     527
     528    psProjectionClass class = PS_PROJ_NONE;
     529    switch (projection->type) {
     530    case PS_PROJ_LIN:
     531    case PS_PROJ_PLY:
     532    case PS_PROJ_WRP:
     533        class = PS_PROJ_CARTESIAN;
     534        break;
     535    case PS_PROJ_TAN:
     536    case PS_PROJ_DIS:
     537    case PS_PROJ_SIN:
     538    case PS_PROJ_STG:
     539    case PS_PROJ_ZEA:
     540    case PS_PROJ_ZPL:
     541        class = PS_PROJ_ZENITHAL;
     542        break;
     543    case PS_PROJ_AIT:
     544    case PS_PROJ_PAR:
     545    case PS_PROJ_GLS:
     546    case PS_PROJ_CAR:
     547    case PS_PROJ_MER:
     548        class = PS_PROJ_PSEUDOCYL;
     549        break;
     550
     551    default:
     552        psAbort ("psCoord", "invalid projection");
     553        break;
     554    }
    455555
    456556    // Allocate return sphere structure
    … …  
    466566    }
    467567
    468     // Remove plate scales
    469     psF64  x = coord->x*projection->Xs;
    470     psF64  y = coord->y*projection->Ys;
    471     psF64  R = sqrt(x*x + y*y);
    472 
    473     bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN);
    474 
    475568    // Perform inverse projection
    476     switch (projection->type) {
    477     case PS_PROJ_TAN:
    478         // Gnonomic deprojection
    479         // the MHPCC version here was fine
    480         # if 0
    481 
     569    switch (class) {
     570    case PS_PROJ_CARTESIAN:
     571        out->r = coord->x*projection->Xs + projection->R;
     572        out->d = coord->y*projection->Ys + projection->D;
     573        break;
     574
     575    case PS_PROJ_ZENITHAL:
     576        // Remove plate scales
     577        x = coord->x*projection->Xs;
     578        y = coord->y*projection->Ys;
     579        R = sqrt(x*x + y*y);
    482580        sinPhi   = (R == 0) ? 0.0 : +x / R;
    483581        cosPhi   = (R == 0) ? 1.0 : -y / R;
    484         rho      = sqrt (1 + R*R);
    485         sinTheta = 1 / rho;
    486         cosTheta = R / rho;
    487         # else
    488 
     582
     583        switch (projection->type)
     584        {
     585        case PS_PROJ_TAN:
     586        case PS_PROJ_DIS:
     587            // Gnonomic deprojection
    489588            rho      = sqrt (1 + R*R);
    490         sinTheta = 1 / rho;
    491         cosTheta = R / rho;
    492         sinPhi   = (R == 0) ? 0.0 : +x / R;
    493         cosPhi   = (R == 0) ? 1.0 : -y / R;
    494         # endif
    495 
    496         break;
    497     case PS_PROJ_SIN:
    498         // Orhtographic deprojection
    499         sinTheta = sqrt (1 - R*R);
    500         cosTheta = R;
    501         sinPhi   = (R == 0) ? 0.0 : +x / R;
    502         cosPhi   = (R == 0) ? 1.0 : -y / R;
    503         break;
    504     case PS_PROJ_AIT:
    505         // Hammer-Aitoff deprojection
    506         // XXX EAM : need range check on z^2 : must be > 0
    507         // XXX EAM : old code, ADD, and elixir code are discrepant re x/4, y/2
    508         rho = sqrt(1.0 - PS_SQR(x/4.0) - PS_SQR(y/2.0));
    509         phi = 2.0*atan2((2.0*rho*rho-1.0), x*rho);
    510         theta = asin(y*rho);
    511         break;
    512     case PS_PROJ_PAR:
    513         // Parabolic deprojection
    514         rho = y/M_PI;
    515         phi = x/(1.0 - 4.0*rho*rho);
    516         theta = 3.0*asin(rho);
    517         break;
    518     default:
    519         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    520                 _("The projection type, %d, is unknown."),
    521                 projection->type);
    522         psFree(out);
    523         return NULL;
    524     }
    525 
    526     if (zenithal) {
     589            sinTheta = 1 / rho;
     590            cosTheta = R / rho;
     591            break;
     592        case PS_PROJ_SIN:
     593            // Orhtographic deprojection
     594            sinTheta = sqrt (1 - R*R);
     595            cosTheta = R;
     596            break;
     597        case PS_PROJ_STG:
     598            psAbort ("coords", "STG not defined");
     599            break;
     600        case PS_PROJ_ZEA:
     601        case PS_PROJ_ZPL:
     602            if (R > 2)
     603                return NULL;
     604            sinTheta = 1 - 0.5*PS_SQR(R);
     605            cosTheta = sqrt (1 - PS_SQR(sinTheta));
     606            break;
     607        default:
     608            psAbort ("psCoord", "invalid projection");
     609            break;
     610        }
     611
    527612        psF64 sinDp = sin(projection->D);
    528613        psF64 cosDp = cos(projection->D);
    … …  
    531616        // psLib versions:
    532617        # if ELIXIR_CODE
    533         // XXX the elixir version : does the ADD have a sign error
    534         psF64 delta = asin(sinTheta*sinDp - cosTheta*cosPhi*cosDp);
    535         psF64 sinAlpha = cosTheta*sinPhi;
    536         psF64 cosAlpha = sinTheta*cosDp + cosTheta*cosPhi*sinDp;
     618        // XXX the elixir version : does the ADD have a sign error?
     619        psF64 delta    = asin(sinTheta*sinDp - cosTheta*cosPhi*cosDp);
     620        psF64 sinAlpha = +cosTheta*sinPhi;
     621        psF64 cosAlpha = +cosTheta*cosPhi*sinDp + sinTheta*cosDp;
    537622        # else
    538623
    539             psF64 delta = asin(sinTheta*sinDp + cosTheta*cosDp*cosPhi);
     624            psF64 delta    = asin(sinTheta*sinDp + cosTheta*cosPhi*cosDp);
    540625        psF64 sinAlpha = -cosTheta*sinPhi;
    541626        psF64 cosAlpha = -cosTheta*cosPhi*sinDp + sinTheta*cosDp;
    … …  
    544629        out->d = delta;
    545630        out->r = atan2(sinAlpha, cosAlpha) + projection->R;
    546     } else {
     631        break;
     632
     633    case PS_PROJ_PSEUDOCYL:
     634        switch (projection->type)
     635        {
     636        case PS_PROJ_AIT:
     637            // Hammer-Aitoff deprojection
     638            // XXX EAM : need range check on z^2 : must be > 0
     639            // XXX EAM : old code, ADD, and elixir code are discrepant re x/4, y/2
     640            rho2 = 1.0 - PS_SQR(0.25*x) - PS_SQR(0.5*y);
     641            if (rho2 < 0)
     642                return (NULL);
     643            rho = sqrt(rho2);
     644            phi = 2.0*atan2(0.5*x*rho, 2.0*rho2 - 1.0);
     645            theta = asin(y*rho);
     646            break;
     647        case PS_PROJ_PAR:
     648            // Parabolic deprojection
     649            rho = y/M_PI;
     650            phi = x/(1.0 - 4.0*rho*rho);
     651            theta = 3.0*asin(rho);
     652            break;
     653        case PS_PROJ_GLS:
     654            phi = x/cos(y);
     655            theta = y;
     656            break;
     657        default:
     658            psAbort ("psCoord", "invalid projection");
     659            break;
     660        }
    547661        out->r = phi   + projection->R;
    548662        out->d = theta + projection->D;
     663        break;
     664    default:
     665        psAbort ("psCoord", "invalid projection");
     666        break;
    549667    }
    550668
    … …  
    10701188    return(sphere);
    10711189}
     1190
     1191char * psProjectTypeToString (psProjectionType type, char *prefix)
     1192{
     1193
     1194    char *name = NULL;
     1195
     1196    switch (type) {
     1197    case PS_PROJ_LIN:
     1198        psStringAppend (&name, "%s-LIN", prefix);
     1199        return name;
     1200    case PS_PROJ_PLY:
     1201        psStringAppend (&name, "%s-PLY", prefix);
     1202        return name;
     1203    case PS_PROJ_WRP:
     1204        psStringAppend (&name, "%s-WRP", prefix);
     1205        return name;
     1206    case PS_PROJ_TAN:
     1207        psStringAppend (&name, "%s-TAN", prefix);
     1208        return name;
     1209    case PS_PROJ_DIS:
     1210        psStringAppend (&name, "%s-DIS", prefix);
     1211        return name;
     1212    case PS_PROJ_SIN:
     1213        psStringAppend (&name, "%s-SIN", prefix);
     1214        return name;
     1215    case PS_PROJ_STG:
     1216        psStringAppend (&name, "%s-STG", prefix);
     1217        return name;
     1218    case PS_PROJ_AIT:
     1219        psStringAppend (&name, "%s-AIT", prefix);
     1220        return name;
     1221    case PS_PROJ_PAR:
     1222        psStringAppend (&name, "%s-PAR", prefix);
     1223        return name;
     1224    case PS_PROJ_GLS:
     1225        psStringAppend (&name, "%s-GLS", prefix);
     1226        return name;
     1227    case PS_PROJ_CAR:
     1228        psStringAppend (&name, "%s-CAR", prefix);
     1229        return name;
     1230    case PS_PROJ_MER:
     1231        psStringAppend (&name, "%s-MER", prefix);
     1232        return name;
     1233
     1234    default:
     1235        psLogMsg ("psLib.astrom", 2, "warning: unknown projection type %d\n", type);
     1236        return NULL;
     1237    }
     1238    return NULL;
     1239}
     1240
     1241psProjectionType psProjectTypeFromString (char *name)
     1242{
     1243
     1244    psProjectionType type = PS_PROJ_NTYPE;
     1245    if (!strcmp (&name[4], "-LIN"))
     1246        type = PS_PROJ_LIN;
     1247    if (!strcmp (&name[4], "-PLY"))
     1248        type = PS_PROJ_PLY;
     1249    if (!strcmp (&name[4], "-WRP"))
     1250        type = PS_PROJ_WRP;
     1251    if (!strcmp (&name[4], "-TAN"))
     1252        type = PS_PROJ_TAN;
     1253    if (!strcmp (&name[4], "-DIS"))
     1254        type = PS_PROJ_DIS;
     1255    if (!strcmp (&name[4], "-SIN"))
     1256        type = PS_PROJ_SIN;
     1257    if (!strcmp (&name[4], "-STG"))
     1258        type = PS_PROJ_STG;
     1259    if (!strcmp (&name[4], "-AIT"))
     1260        type = PS_PROJ_AIT;
     1261    if (!strcmp (&name[4], "-PAR"))
     1262        type = PS_PROJ_PAR;
     1263    if (!strcmp (&name[4], "-GLS"))
     1264        type = PS_PROJ_GLS;
     1265    if (!strcmp (&name[4], "-CAR"))
     1266        type = PS_PROJ_CAR;
     1267    if (!strcmp (&name[4], "-MER"))
     1268        type = PS_PROJ_MER;
     1269    return type;
     1270}

  • trunk/psLib/src/astro/psCoord.h

    r10609 r10823  
    1111 *  @author GLG, MHPCC
    1212 *
    13  *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2006-12-10 17:29:51 $
     13 *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2006-12-22 21:19:47 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    122122 */
    123123typedef enum {
     124    PS_PROJ_LIN,                       ///< Linear projection
     125    PS_PROJ_PLY,                       ///< Linear polynomial projection
     126    PS_PROJ_WRP,                       ///< Linear polynomial projection
    124127    PS_PROJ_TAN,                       ///< Tangent projection
     128    PS_PROJ_DIS,                       ///< Sine projection
    125129    PS_PROJ_SIN,                       ///< Sine projection
     130    PS_PROJ_STG,                       ///< Sine projection
     131    PS_PROJ_ZEA,                       ///< Sine projection
     132    PS_PROJ_ZPL,                       ///< Sine projection
    126133    PS_PROJ_AIT,                       ///< Aitoff projection
    127134    PS_PROJ_PAR,                       ///< Par projection
    128     //    PS_PROJ_GLS,                ///< GLS projection
    129     //    PS_PROJ_CAR,                ///< CAR projection
    130     //    PS_PROJ_MER,                ///< MER projection
     135    PS_PROJ_GLS,                ///< GLS projection
     136    PS_PROJ_CAR,                ///< CAR projection
     137    PS_PROJ_MER,                ///< MER projection
    131138    PS_PROJ_NTYPE                      ///< Number of types; must be last.
    132139} psProjectionType;
    … …  
    438445);
    439446
     447char * psProjectTypeToString (psProjectionType type, char *prefix);
     448psProjectionType psProjectTypeFromString (char *name);
     449
     450
    440451#define PS_PRINT_PLANE_TRANSFORM(NAME) \
    441452{ \
Note: See TracChangeset for help on using the changeset viewer.