Application: gvSIG desktop

/*

** libproj -- library of cartographic projections

**

** Copyright (c) 2003   Gerald I. Evenden

*/

static const char

LIBPROJ_ID[] = "$Id: pj_gauss.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 PJ_LIB__

#include <projects.h>

#define MAX_ITER 20

struct GAUSS {

        double C;

        double K;

        double e;

        double ratexp;

};

#define EN ((struct GAUSS *)en)

#define DEL_TOL 1e-14

        static double

srat(double esinp, double exp) {

        return(pow((1.-esinp)/(1.+esinp), exp));

}

        void *

pj_gauss_ini(double e, double phi0, double *chi, double *rc) {

        double sphi, cphi, es;

        struct GAUSS *en;

        if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL)

                return (NULL);

        es = e * e;

        EN->e = e;

        sphi = sin(phi0);

        cphi = cos(phi0);  cphi *= cphi;

        *rc = sqrt(1. - es) / (1. - es * sphi * sphi);

        EN->C = sqrt(1. + es * cphi * cphi / (1. - es));

        *chi = asin(sphi / EN->C);

        EN->ratexp = 0.5 * EN->C * e;

        EN->K = tan(.5 * *chi + FORTPI) / (

                pow(tan(.5 * phi0 + FORTPI), EN->C) *

                srat(EN->e * sphi, EN->ratexp)  );

        return ((void *)en);

}

        LP

pj_gauss(LP elp, const void *en) {

        LP slp;

        slp.phi = 2. * atan( EN->K *

                pow(tan(.5 * elp.phi + FORTPI), EN->C) *

                srat(EN->e * sin(elp.phi), EN->ratexp) ) - HALFPI;

        slp.lam = EN->C * (elp.lam);

        return(slp);

}

        LP

pj_inv_gauss(LP slp, const void *en) {

        LP elp;

        double num;

        int i;

        elp.lam = slp.lam / EN->C;

        num = pow(tan(.5 * slp.phi + FORTPI)/EN->K, 1./EN->C);

        for (i = MAX_ITER; i; --i) {

                elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e))

                        - HALFPI;

                if (fabs(elp.phi - slp.phi) < DEL_TOL) break;

                        slp.phi = elp.phi;

        }        

        /* convergence failed */

        if (!i)

                pj_errno = -17;

        return (elp);

}

/* Revision Log:

** $Log: pj_gauss.c,v $

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

** New

**

** Revision 2.2  2004/03/15 16:07:42  gie

** removed es from init structure

**

** Revision 2.1  2003/03/28 01:44:30  gie

** Initial

**

*/