Application: gvSIG desktop

/*

** libproj -- library of cartographic projections

**

** Copyright (c) 2004   Gerald I. Evenden

*/

static const char

LIBPROJ_ID[] = "$Id: PJ_geos.c,v 1.1 2004/10/20 17:04:00 fwarmerdam Exp $";

/*

** Permission is hereby granted, free of charge, to any person obtaining

** a copy of this software and associated documentation files (the

** "Software"), to deal in the Software without restriction, including

** without limitation the rights to use, copy, modify, merge, publish,

** distribute, sublicense, and/or sell copies of the Software, and to

** permit persons to whom the Software is furnished to do so, subject to

** the following conditions:

**

** The above copyright notice and this permission notice shall be

** included in all copies or substantial portions of the Software.

**

** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,

** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF

** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.

** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY

** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,

** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE

** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

*/

#define PROJ_PARMS__ \

        double        h; \

        double  radius_p; \

        double  radius_p2; \

        double  radius_p_inv2; \

        double  radius_g; \

        double  radius_g_1; \

        double  C;

#define PJ_LIB__

#include        <projects.h>

PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th=";

FORWARD(s_forward); /* spheroid */

        double Vx, Vy, Vz, tmp;

/* Calculation of the three components of the vector from satellite to

** position on earth surface (lon,lat).*/

        tmp = cos(lp.phi);

        Vx = cos (lp.lam) * tmp;

        Vy = sin (lp.lam) * tmp;

        Vz = sin (lp.phi);

/* Check visibility.*/

        if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR;

/* Calculation based on view angles from satellite.*/

        tmp = P->radius_g - Vx;

        xy.x = P->radius_g_1 * atan(Vy / tmp);

        xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));

        return (xy);

}

FORWARD(e_forward); /* ellipsoid */

        double r, Vx, Vy, Vz, tmp;

/* Calculation of geocentric latitude. */

        lp.phi = atan (P->radius_p2 * tan (lp.phi));

/* Calculation of the three components of the vector from satellite to

** position on earth surface (lon,lat).*/

        r = (P->radius_p) / hypot(P->radius_p * cos (lp.phi), sin (lp.phi));

        Vx = r * cos (lp.lam) * cos (lp.phi);

        Vy = r * sin (lp.lam) * cos (lp.phi);

        Vz = r * sin (lp.phi);

/* Check visibility. */

        if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * P->radius_p_inv2) < 0.)

                F_ERROR;

/* Calculation based on view angles from satellite. */

        tmp = P->radius_g - Vx;

        xy.x = P->radius_g_1 * atan (Vy / tmp);

        xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));

        return (xy);

}

INVERSE(s_inverse); /* spheroid */

        double Vx, Vy, Vz, a, b, c, det, k;

/* Setting three components of vector from satellite to position.*/

        Vx = -1.0;

        Vy = tan (xy.x / (P->radius_g - 1.0));

        Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);

/* Calculation of terms in cubic equation and determinant.*/

        a   = Vy * Vy + Vz * Vz + Vx * Vx;

        b   = 2 * P->radius_g * Vx;

        if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;

/* Calculation of three components of vector from satellite to position.*/

        k  = (-b - sqrt(det)) / (2 * a);

        Vx = P->radius_g + k * Vx;

        Vy *= k;

        Vz *= k;

/* Calculation of longitude and latitude.*/

        lp.lam = atan2 (Vy, Vx);

        lp.phi = atan (Vz * cos (lp.lam) / Vx);

        return (lp);

}

INVERSE(e_inverse); /* ellipsoid */

        double Vx, Vy, Vz, a, b, c, det, k;

/* Setting three components of vector from satellite to position.*/

        Vx = -1.0;

        Vy = tan (xy.x / P->radius_g_1);

        Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);

/* Calculation of terms in cubic equation and determinant.*/

        a = Vz / P->radius_p;

        a   = Vy * Vy + a * a + Vx * Vx;

        b   = 2 * P->radius_g * Vx;

        if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;

/* Calculation of three components of vector from satellite to position.*/

        k  = (-b - sqrt(det)) / (2. * a);

        Vx = P->radius_g + k * Vx;

        Vy *= k;

        Vz *= k;

/* Calculation of longitude and latitude.*/

        lp.lam  = atan2 (Vy, Vx);

        lp.phi = atan (Vz * cos (lp.lam) / Vx);

        lp.phi = atan (P->radius_p_inv2 * tan (lp.phi));

        return (lp);

}

FREEUP; if (P) free(P); }

ENTRY0(geos)

        if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);

        if (P->phi0) E_ERROR(-46);

        P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);

        P->C  = P->radius_g * P->radius_g - 1.0;

        if (P->es) {

                P->radius_p      = sqrt (P->one_es);

                P->radius_p2     = P->one_es;

                P->radius_p_inv2 = P->rone_es;

                P->inv = e_inverse;

                P->fwd = e_forward;

        } else {

                P->radius_p = P->radius_p2 = P->radius_p_inv2 = 1.0;

                P->inv = s_inverse;

                P->fwd = s_forward;

        }

ENDENTRY(P)

/*

** $Log: PJ_geos.c,v $

** Revision 1.1  2004/10/20 17:04:00  fwarmerdam

** New

**

** Revision 1.2  2004/07/14 18:08:57  gie

** corrected P->phi_0 to P->phi0

**

** Revision 1.1  2004/07/12 17:58:25  gie

** Initial revision

**

*/