IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3450


Ignore:
Timestamp:
Mar 18, 2005, 11:34:43 AM (21 years ago)
Author:
evanalst
Message:

Update psProject and psDeproject functions to match ADD-09.

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r3359 r3450  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-02 03:39:58 $
     12*  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-18 21:34:43 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    399399}
    400400
    401 // XXX: Must test psProjectionAlloc() and projectionFree().
    402 // XXX: Must rewrite code and tests to use these functions.
    403401psProjection* psProjectionAlloc(
    404402    psF64 R,
    … …  
    419417}
    420418
    421 
    422 /******************************************************************************
    423 XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    424 XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    425 XXX: Must apply scaling at the end.
    426  *****************************************************************************/
    427419psPlane* psProject(const psSphere* coord,
    428420                   const psProjection* projection)
    … …  
    431423    PS_PTR_CHECK_NULL(projection, NULL);
    432424
    433     psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
     425    psF64   theta = 0.0;
     426    psF64   phi   = 0.0;
     427
     428    // Allocate return value
     429    psPlane* out = psPlaneAlloc();
     430
     431    // Convert to projection spherical coordinate system
     432    theta = asin( sin(coord->d)*sin(projection->D) +
     433                  cos(coord->d)*cos(projection->D)*cos(coord->r-projection->R));
     434    phi = atan2( -1.0*cos(coord->d)*sin(coord->r-projection->R),
     435                 sin(coord->d)*cos(projection->D) - cos(coord->d)*sin(projection->D)*cos(coord->r-projection->R) );
     436
     437    // Perform the specified projection
     438    // Gnomonic projection
    434439    if (projection->type == PS_PROJ_TAN) {
    435 
    436         /*
    437                 // ********************************************
    438                 // From the ADD
    439                 // ********************************************
    440                 // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
    441                 // return(NULL);
    442                 psF32 sinTheta;
    443                 psF32 cosThetaCosPhi;
    444                 psF32 cosThetaSinPhi;
    445          
    446         sinTheta = (sin(coord->d) * sin(projection->D)) +
    447              (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R)));
    448         cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) -
    449                    (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R)));
    450         cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R);
    451          
    452                 tmp->x =  -cosThetaSinPhi / sinTheta;
    453                 tmp->y = cosThetaCosPhi / sinTheta;
    454          
    455         */
    456         // ********************************************
    457         // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
    458         // delta_0 and phi_1 are the projection centers, not the point being projected.
    459         // (delta_0, phi_1) == (projection->R, projection->D)
    460         // ********************************************
    461         psF32 cosC = (sin(projection->D) * sin(coord->d)) +
    462                      (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R));
    463         tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC;
    464         tmp->y = ((cos(projection->D) * sin(coord->d)) -
    465                   (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) /
    466                  cosC;
    467 
     440        out->x = (cos(theta)*sin(phi))/sin(theta);
     441        out->y = (-1.0*cos(theta)*cos(phi))/sin(theta);
     442        // Othrographic projection
    468443    } else if (projection->type == PS_PROJ_SIN) {
    469         // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
    470         // return(NULL);
    471 
    472         tmp->x = cos(coord->d) * sin(coord->d);
    473         tmp->y = -cos(coord->d) * cos(coord->d);
    474 
    475     } else if (projection->type == PS_PROJ_CAR) {
    476         tmp->x = coord->d;
    477         tmp->y = coord->r;
    478 
    479     } else if (projection->type == PS_PROJ_MER) {
    480         tmp->x = coord->d;
    481         tmp->y = log(tan(DEG_TO_RAD(45.0) + (0.5 * coord->r)));
    482     } else if (projection->type == PS_PROJ_AIT) {
    483         psF64 tmpF64  = PS_SQRT_F32(0.5 * (1.0 + (cos(coord->r) * cos(0.5 * coord->d))));
    484         tmpF64 = 1.0 / tmpF64;
    485         tmp->x = 2.0 * tmpF64 * cos(coord->r) * sin(0.5 * coord->d);
    486         tmp->y = sin(coord->r) * tmpF64;
    487 
    488     } else if (projection->type == PS_PROJ_PAR) {
    489         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    490                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    491                 "PS_PROJ_PAR");
    492 
    493     } else if (projection->type == PS_PROJ_GLS) {
    494         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    495                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    496                 "PS_PROJ_GLS");
     444        out->x = cos(theta)*sin(phi);
     445        out->y = -1.0*cos(theta)*cos(phi);
     446        // Hammer-Aitoff projection
     447    } else if ( projection->type == PS_PROJ_AIT) {
     448        psF64 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0)));
     449        out->x = 2.0*zeta*cos(theta)*sin(phi/2.0);
     450        out->y = zeta*sin(theta);
     451        // Parabolic projection
     452    } else if ( projection->type == PS_PROJ_PAR) {
     453        out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0);
     454        out->y = PS_PI*sin(theta/3.0);
    497455    } else {
    498456        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    499457                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
    500458                projection->type);
    501     }
    502 
    503     return (tmp);
    504 }
    505 
    506 /******************************************************************************
    507 XXX: Waiting for algorithms.
    508 XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    509 XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    510 XXX: Must apply scaling at the beginning.
    511  *****************************************************************************/
     459        psFree(out);
     460        return NULL;
     461    }
     462
     463    // Apply plate scales
     464    out->x *= projection->Xs;
     465    out->y *= projection->Ys;
     466
     467    // Return output
     468    return out;
     469}
     470
    512471psSphere* psDeproject(const psPlane* coord,
    513472                      const psProjection* projection)
    … …  
    515474    PS_PTR_CHECK_NULL(coord, NULL);
    516475    PS_PTR_CHECK_NULL(projection, NULL);
    517     //    psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work.");
    518     //    return(NULL);
    519     float R = 0.0;
    520     float chu = 0.0;
    521     float chu1 = 0.0;
    522     float chu2 = 0.0;
    523     psSphere* tmp = (psSphere* ) psAlloc(sizeof(psSphere));
    524 
    525     if (projection->type == PS_PROJ_TAN) {
    526         /*
    527                 // ********************************************
    528                 // From the ADD
    529                 // ********************************************
    530                 psF32 phi;
    531                 psF32 theta;
    532                 psF32 sinDelta;
    533                 psF32 cosDeltaCosAlphaMinusAlphaP;
    534                 psF32 cosDeltaSinAlphaMinusAlphaP;
    535          
    536                 phi = atan2(-coord->y, coord->x);
    537                 R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    538                 theta = atan(1.0 / R);
    539          
    540                 sinDelta = (sin(theta) * sin(projection->D)) +
    541                            (cos(theta) * cos(projection->D) * cos(phi));
    542                 cosDeltaCosAlphaMinusAlphaP =
    543                            (sin(theta) * cos(projection->D)) -
    544                            (cos(theta) * sin(projection->D) * cos(phi));
    545                 cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi);
    546                 tmp->d = asin(sinDelta);
    547                 tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R;
    548         */
    549         // ********************************************
    550         // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
    551         // delta_0 and phi_1 are the projection centers, not the point being projected.
    552         // (delta_0, phi_1) == (projection->R, projection->D)
    553         // XXX: figure out how to use the two-argument atan2() here.
    554         // ********************************************
    555         psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    556         psF32 C = atan(row);
    557         tmp->d = asin((cos(C) * sin(projection->D)) +
    558                       ((coord->y * sin(C) * cos(projection->D)) / row));
    559         psF32 tmpAtan = atan((coord->x * sin(C)) /
    560                              ((row * cos(projection->D) * cos(C)) -
    561                               (coord->y * sin(projection->D) * sin(C))));
    562         tmp->r = projection->R + tmpAtan;
    563 
    564 
    565     } else if (projection->type == PS_PROJ_SIN) {
    566         R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    567         tmp->d = atan2(-coord->y, coord->x);
    568         tmp->r = acos((R * PS_PI) / 180.0);
    569 
    570     } else if (projection->type == PS_PROJ_CAR) {
    571         tmp->d = coord->x;
    572         tmp->r = coord->y;
    573 
    574     } else if (projection->type == PS_PROJ_MER) {
    575         tmp->d = coord->x;
    576         tmp->r = (2.0 * atan(exp((coord->y * PS_PI / 180.0)))) - 180.0;
    577 
    578     } else if (projection->type == PS_PROJ_AIT) {
    579         chu1 = (coord->x * PS_PI) / 720.0;
    580         chu1 *= chu1;
    581         chu2 = (coord->y * PS_PI) / 360.0;
    582         chu2 *= chu2;
    583         chu = PS_SQRT_F32(1.0 - chu1 - chu2);
    584         tmp->d = 2.0 * atan2((2.0 * chu * chu) - 1.0, (coord->x * chu * PS_PI) / 360.0);
    585         tmp->r = asin((coord->y * chu * PS_PI) / 180.0);
    586 
    587     } else if (projection->type == PS_PROJ_PAR) {
    588         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    589                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    590                 "PS_PROJ_PAR");
    591 
    592     } else if (projection->type == PS_PROJ_GLS) {
    593         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    594                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    595                 "PS_PROJ_GLS");
     476
     477    psF64  theta = 0.0;
     478    psF64  phi   = 0.0;
     479
     480    // Allocate return sphere structure
     481    psSphere* out = psSphereAlloc();
     482
     483    // Remove plate scales
     484    psF64  x = coord->x/projection->Xs;
     485    psF64  y = coord->y/projection->Ys;
     486
     487    // Perform inverse projection
     488    // Gnonomic deprojection
     489    if ( projection->type == PS_PROJ_TAN) {
     490        phi = atan(-1.0*x/y);
     491        theta = atan(1.0/sqrt(x*x+y*y));
     492        // Orhtographic deprojection
     493    } else if ( projection->type == PS_PROJ_SIN) {
     494        phi = atan((-1.0*x)/y);
     495        theta = atan( sqrt(1.0-(x*x+y*y)) / sqrt(x*x+y*y));
     496        // Hammer-Aitoff deprojection
     497    } else if ( projection->type == PS_PROJ_AIT) {
     498        psF64 z = sqrt(1.0 - ((x/4.0)*(x/4.0)) - ((y/2.0)*(y/2.0)));
     499        phi = 2.0*atan((z*x) / (2.0*(2.0*z*z-1.0)) );
     500        theta = asin(y*z);
     501        // Parabolic deprojection
     502    } else if ( projection->type == PS_PROJ_PAR) {
     503        psF64 rho = y/PS_PI;
     504        phi = x/(1.0 - 4.0*rho*rho);
     505        theta = 3.0*asin(rho);
     506        // Invalid deprojection type
    596507    } else {
    597508        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    598509                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
    599510                projection->type);
    600     }
    601 
    602     return (tmp);
     511        psFree(out);
     512        return NULL;
     513    }
     514
     515    // Convert from projection spherical coordinates
     516    out->d = asin( sin(theta)*sin(projection->D) +
     517                   cos(theta)*cos(projection->D)*cos(phi) );
     518    out->r = projection->R + atan2( -1.0*cos(theta)*sin(phi),
     519                                    sin(theta)*cos(projection->D) -
     520                                    cos(theta)*sin(projection->D)*cos(phi) );
     521
     522    // Return sphere coordinate
     523    return out;
    603524}
    604525

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

    r3334 r3450  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-02-25 22:29:19 $
     12*  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-18 21:34:43 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    124124    PS_PROJ_AIT,                ///< Aitoff projection
    125125    PS_PROJ_PAR,                ///< Par projection
    126     PS_PROJ_GLS,                ///< GLS projection
    127     PS_PROJ_CAR,                ///< CAR projection
    128     PS_PROJ_MER,                ///< MER projection
     126    //    PS_PROJ_GLS,                ///< GLS projection
     127    //    PS_PROJ_CAR,                ///< CAR projection
     128    //    PS_PROJ_MER,                ///< MER projection
    129129    PS_PROJ_NTYPE               ///< Number of types; must be last.
    130130} psProjectionType;

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

    r3359 r3450  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-02 03:39:58 $
     12*  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-18 21:34:43 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    399399}
    400400
    401 // XXX: Must test psProjectionAlloc() and projectionFree().
    402 // XXX: Must rewrite code and tests to use these functions.
    403401psProjection* psProjectionAlloc(
    404402    psF64 R,
    … …  
    419417}
    420418
    421 
    422 /******************************************************************************
    423 XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    424 XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    425 XXX: Must apply scaling at the end.
    426  *****************************************************************************/
    427419psPlane* psProject(const psSphere* coord,
    428420                   const psProjection* projection)
    … …  
    431423    PS_PTR_CHECK_NULL(projection, NULL);
    432424
    433     psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
     425    psF64   theta = 0.0;
     426    psF64   phi   = 0.0;
     427
     428    // Allocate return value
     429    psPlane* out = psPlaneAlloc();
     430
     431    // Convert to projection spherical coordinate system
     432    theta = asin( sin(coord->d)*sin(projection->D) +
     433                  cos(coord->d)*cos(projection->D)*cos(coord->r-projection->R));
     434    phi = atan2( -1.0*cos(coord->d)*sin(coord->r-projection->R),
     435                 sin(coord->d)*cos(projection->D) - cos(coord->d)*sin(projection->D)*cos(coord->r-projection->R) );
     436
     437    // Perform the specified projection
     438    // Gnomonic projection
    434439    if (projection->type == PS_PROJ_TAN) {
    435 
    436         /*
    437                 // ********************************************
    438                 // From the ADD
    439                 // ********************************************
    440                 // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
    441                 // return(NULL);
    442                 psF32 sinTheta;
    443                 psF32 cosThetaCosPhi;
    444                 psF32 cosThetaSinPhi;
    445          
    446         sinTheta = (sin(coord->d) * sin(projection->D)) +
    447              (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R)));
    448         cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) -
    449                    (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R)));
    450         cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R);
    451          
    452                 tmp->x =  -cosThetaSinPhi / sinTheta;
    453                 tmp->y = cosThetaCosPhi / sinTheta;
    454          
    455         */
    456         // ********************************************
    457         // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
    458         // delta_0 and phi_1 are the projection centers, not the point being projected.
    459         // (delta_0, phi_1) == (projection->R, projection->D)
    460         // ********************************************
    461         psF32 cosC = (sin(projection->D) * sin(coord->d)) +
    462                      (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R));
    463         tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC;
    464         tmp->y = ((cos(projection->D) * sin(coord->d)) -
    465                   (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) /
    466                  cosC;
    467 
     440        out->x = (cos(theta)*sin(phi))/sin(theta);
     441        out->y = (-1.0*cos(theta)*cos(phi))/sin(theta);
     442        // Othrographic projection
    468443    } else if (projection->type == PS_PROJ_SIN) {
    469         // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
    470         // return(NULL);
    471 
    472         tmp->x = cos(coord->d) * sin(coord->d);
    473         tmp->y = -cos(coord->d) * cos(coord->d);
    474 
    475     } else if (projection->type == PS_PROJ_CAR) {
    476         tmp->x = coord->d;
    477         tmp->y = coord->r;
    478 
    479     } else if (projection->type == PS_PROJ_MER) {
    480         tmp->x = coord->d;
    481         tmp->y = log(tan(DEG_TO_RAD(45.0) + (0.5 * coord->r)));
    482     } else if (projection->type == PS_PROJ_AIT) {
    483         psF64 tmpF64  = PS_SQRT_F32(0.5 * (1.0 + (cos(coord->r) * cos(0.5 * coord->d))));
    484         tmpF64 = 1.0 / tmpF64;
    485         tmp->x = 2.0 * tmpF64 * cos(coord->r) * sin(0.5 * coord->d);
    486         tmp->y = sin(coord->r) * tmpF64;
    487 
    488     } else if (projection->type == PS_PROJ_PAR) {
    489         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    490                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    491                 "PS_PROJ_PAR");
    492 
    493     } else if (projection->type == PS_PROJ_GLS) {
    494         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    495                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    496                 "PS_PROJ_GLS");
     444        out->x = cos(theta)*sin(phi);
     445        out->y = -1.0*cos(theta)*cos(phi);
     446        // Hammer-Aitoff projection
     447    } else if ( projection->type == PS_PROJ_AIT) {
     448        psF64 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0)));
     449        out->x = 2.0*zeta*cos(theta)*sin(phi/2.0);
     450        out->y = zeta*sin(theta);
     451        // Parabolic projection
     452    } else if ( projection->type == PS_PROJ_PAR) {
     453        out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0);
     454        out->y = PS_PI*sin(theta/3.0);
    497455    } else {
    498456        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    499457                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
    500458                projection->type);
    501     }
    502 
    503     return (tmp);
    504 }
    505 
    506 /******************************************************************************
    507 XXX: Waiting for algorithms.
    508 XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    509 XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    510 XXX: Must apply scaling at the beginning.
    511  *****************************************************************************/
     459        psFree(out);
     460        return NULL;
     461    }
     462
     463    // Apply plate scales
     464    out->x *= projection->Xs;
     465    out->y *= projection->Ys;
     466
     467    // Return output
     468    return out;
     469}
     470
    512471psSphere* psDeproject(const psPlane* coord,
    513472                      const psProjection* projection)
    … …  
    515474    PS_PTR_CHECK_NULL(coord, NULL);
    516475    PS_PTR_CHECK_NULL(projection, NULL);
    517     //    psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work.");
    518     //    return(NULL);
    519     float R = 0.0;
    520     float chu = 0.0;
    521     float chu1 = 0.0;
    522     float chu2 = 0.0;
    523     psSphere* tmp = (psSphere* ) psAlloc(sizeof(psSphere));
    524 
    525     if (projection->type == PS_PROJ_TAN) {
    526         /*
    527                 // ********************************************
    528                 // From the ADD
    529                 // ********************************************
    530                 psF32 phi;
    531                 psF32 theta;
    532                 psF32 sinDelta;
    533                 psF32 cosDeltaCosAlphaMinusAlphaP;
    534                 psF32 cosDeltaSinAlphaMinusAlphaP;
    535          
    536                 phi = atan2(-coord->y, coord->x);
    537                 R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    538                 theta = atan(1.0 / R);
    539          
    540                 sinDelta = (sin(theta) * sin(projection->D)) +
    541                            (cos(theta) * cos(projection->D) * cos(phi));
    542                 cosDeltaCosAlphaMinusAlphaP =
    543                            (sin(theta) * cos(projection->D)) -
    544                            (cos(theta) * sin(projection->D) * cos(phi));
    545                 cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi);
    546                 tmp->d = asin(sinDelta);
    547                 tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R;
    548         */
    549         // ********************************************
    550         // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
    551         // delta_0 and phi_1 are the projection centers, not the point being projected.
    552         // (delta_0, phi_1) == (projection->R, projection->D)
    553         // XXX: figure out how to use the two-argument atan2() here.
    554         // ********************************************
    555         psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    556         psF32 C = atan(row);
    557         tmp->d = asin((cos(C) * sin(projection->D)) +
    558                       ((coord->y * sin(C) * cos(projection->D)) / row));
    559         psF32 tmpAtan = atan((coord->x * sin(C)) /
    560                              ((row * cos(projection->D) * cos(C)) -
    561                               (coord->y * sin(projection->D) * sin(C))));
    562         tmp->r = projection->R + tmpAtan;
    563 
    564 
    565     } else if (projection->type == PS_PROJ_SIN) {
    566         R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    567         tmp->d = atan2(-coord->y, coord->x);
    568         tmp->r = acos((R * PS_PI) / 180.0);
    569 
    570     } else if (projection->type == PS_PROJ_CAR) {
    571         tmp->d = coord->x;
    572         tmp->r = coord->y;
    573 
    574     } else if (projection->type == PS_PROJ_MER) {
    575         tmp->d = coord->x;
    576         tmp->r = (2.0 * atan(exp((coord->y * PS_PI / 180.0)))) - 180.0;
    577 
    578     } else if (projection->type == PS_PROJ_AIT) {
    579         chu1 = (coord->x * PS_PI) / 720.0;
    580         chu1 *= chu1;
    581         chu2 = (coord->y * PS_PI) / 360.0;
    582         chu2 *= chu2;
    583         chu = PS_SQRT_F32(1.0 - chu1 - chu2);
    584         tmp->d = 2.0 * atan2((2.0 * chu * chu) - 1.0, (coord->x * chu * PS_PI) / 360.0);
    585         tmp->r = asin((coord->y * chu * PS_PI) / 180.0);
    586 
    587     } else if (projection->type == PS_PROJ_PAR) {
    588         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    589                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    590                 "PS_PROJ_PAR");
    591 
    592     } else if (projection->type == PS_PROJ_GLS) {
    593         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    594                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
    595                 "PS_PROJ_GLS");
     476
     477    psF64  theta = 0.0;
     478    psF64  phi   = 0.0;
     479
     480    // Allocate return sphere structure
     481    psSphere* out = psSphereAlloc();
     482
     483    // Remove plate scales
     484    psF64  x = coord->x/projection->Xs;
     485    psF64  y = coord->y/projection->Ys;
     486
     487    // Perform inverse projection
     488    // Gnonomic deprojection
     489    if ( projection->type == PS_PROJ_TAN) {
     490        phi = atan(-1.0*x/y);
     491        theta = atan(1.0/sqrt(x*x+y*y));
     492        // Orhtographic deprojection
     493    } else if ( projection->type == PS_PROJ_SIN) {
     494        phi = atan((-1.0*x)/y);
     495        theta = atan( sqrt(1.0-(x*x+y*y)) / sqrt(x*x+y*y));
     496        // Hammer-Aitoff deprojection
     497    } else if ( projection->type == PS_PROJ_AIT) {
     498        psF64 z = sqrt(1.0 - ((x/4.0)*(x/4.0)) - ((y/2.0)*(y/2.0)));
     499        phi = 2.0*atan((z*x) / (2.0*(2.0*z*z-1.0)) );
     500        theta = asin(y*z);
     501        // Parabolic deprojection
     502    } else if ( projection->type == PS_PROJ_PAR) {
     503        psF64 rho = y/PS_PI;
     504        phi = x/(1.0 - 4.0*rho*rho);
     505        theta = 3.0*asin(rho);
     506        // Invalid deprojection type
    596507    } else {
    597508        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    598509                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
    599510                projection->type);
    600     }
    601 
    602     return (tmp);
     511        psFree(out);
     512        return NULL;
     513    }
     514
     515    // Convert from projection spherical coordinates
     516    out->d = asin( sin(theta)*sin(projection->D) +
     517                   cos(theta)*cos(projection->D)*cos(phi) );
     518    out->r = projection->R + atan2( -1.0*cos(theta)*sin(phi),
     519                                    sin(theta)*cos(projection->D) -
     520                                    cos(theta)*sin(projection->D)*cos(phi) );
     521
     522    // Return sphere coordinate
     523    return out;
    603524}
    604525

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

    r3334 r3450  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-02-25 22:29:19 $
     12*  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-18 21:34:43 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    124124    PS_PROJ_AIT,                ///< Aitoff projection
    125125    PS_PROJ_PAR,                ///< Par projection
    126     PS_PROJ_GLS,                ///< GLS projection
    127     PS_PROJ_CAR,                ///< CAR projection
    128     PS_PROJ_MER,                ///< MER projection
     126    //    PS_PROJ_GLS,                ///< GLS projection
     127    //    PS_PROJ_CAR,                ///< CAR projection
     128    //    PS_PROJ_MER,                ///< MER projection
    129129    PS_PROJ_NTYPE               ///< Number of types; must be last.
    130130} psProjectionType;
Note: See TracChangeset for help on using the changeset viewer.