#include    <stdio.h>
#include    <stdlib.h>
#include    <math.h>

#define RAD (180.*60./3.14159265359)
#define HPI (3.14159265359/2.)
#define CR  (1.e-9)

/************************************************************************/
/* <<< Coordinate (Projection) Number (NC) >>>                          */
/*----------------------------------------------------------------------*/
/*  For TOKYO Geodetic Datum  in reference to Bessel's ellipsoid        */
/*    0  Japanese Transverse Mercator Proj. (XY)               [0.9999] */
/* 1-60  Universal Transverse Mercator Proj. (UTM) Zone 1-60   [0.9996] */
/*   61  Universal Polar Stereographic Proj. (UPS) North-Pole  [0.994]  */
/*   62  Universal Polar Stereographic Proj. (UPS) South-Pole  [0.994]  */
/*   65  same as UTM but non-standard Central Meridian         [0.9996] */
/*   70  Mercator Projection (Regular)                                  */
/*   71  Lambert Conformal Conic Proj. (1 St.Parallel)                  */
/*   72  Lambert Conformal Conic Proj. (2 St.Parallels)                 */
/*  100  Lambert Azimuthal Equal-Area Projection  (R = Same surf.area)  */
/*  109  Lambert Azimuthal Equal-Area Projection  (R = Semi major axis) */
/*  199  Latitude/Longitude in minutes are regarded as distance (in km) */
/*----------------------------------------------------------------------*/
/*  For World Geodetic Datum  in reference to GWS80 ellipsoid           */
/*  + 200   (Same as above plus 200)                                    */
/************************************************************************/

/******** Bessel's ellipsoid ********/
#define AB  (6377.397155)
#define EB  (0.0816968312)
#define EEB (0.006674372227)
#define E1B (0.001674184800)
/******** GRS80 ellipsoid ***********/
#define AW  (6378.137)
#define EW  (0.081819191)
#define EEW (0.0066943800)
#define E1W (0.0016792204)
/************************************/
/*  Ac  : AB / AW   : (Earth's equatorial radius)                       */
/*  Ec  : EB / EW   : (eccentricity) ;   EEc = Ec*Ec                    */
/*  E1c : E1B / E1W : ((1.-sqrt(1.-EEc)) / (1.+sqrt(1.-EEc)))           */
/*  SRc : SRB / SRW : (radius of sphere with same surface area)         */
/************************************************************************/

#define E4B (EEB*EEB)
#define E6B (EEB*EEB*EEB)
#define E8B (EEB*EEB*EEB*EEB)
#define FRB (1. - EEB/6. - E4B*(17./360.) - E6B*(187./9072.))
#define SRB (AB*FRB)
#define EHB (EB/2.)
#define A2B (AB*1.414213562373)
#define SR2B (SRB*1.414213562373)
#define C0B (AB*(1. -EEB/4. -E4B*(3./64.) -E6B*(5./256.) -E8B*(175./16384.)))
#define C1B (AB*(EEB*(3./4.)+E4B*(27./64.)+E6B*(85./256.)+E8B*(4655./16384.)))
#define C2B (AB*(E4B*(15./32.) + E6B*(275./384.) + E8B*(22225./24576.)))
#define C3B (AB*(E6B*(35./96.) + E8B*(5635./6144.)))
#define F1B (E1B*3. - E1B*E1B*(21./4.) + E1B*E1B*E1B*(31./4.))
#define F2B (E1B*E1B*(21./2.) - E1B*E1B*E1B*(151./3.))
#define F3B (E1B*E1B*E1B*(151./3.))

#define E4W (EEW*EEW)
#define E6W (EEW*EEW*EEW)
#define E8W (EEW*EEW*EEW*EEW)
#define FRW (1. - EEW/6. - E4W*(17./360.) - E6W*(187./9072.))
#define SRW (AW*FRW)
#define EHW (EW/2.)
#define A2W (AW*1.414213562373)
#define SR2W (SRW*1.414213562373)
#define C0W (AW*(1. -EEW/4. -E4W*(3./64.) -E6W*(5./256.) -E8W*(175./16384.)))
#define C1W (AW*(EEW*(3./4.)+E4W*(27./64.)+E6W*(85./256.)+E8W*(4655./16384.)))
#define C2W (AW*(E4W*(15./32.) + E6W*(275./384.) + E8W*(22225./24576.)))
#define C3W (AW*(E6W*(35./96.) + E8W*(5635./6144.)))
#define F1W (E1W*3. - E1W*E1W*(21./4.) + E1W*E1W*E1W*(31./4.))
#define F2W (E1W*E1W*(21./2.) - E1W*E1W*E1W*(151./3.))
#define F3W (E1W*E1W*E1W*(151./3.))
static int nc=0, kw=0;
static double glon, gx, ct, st;
static double Ac, Ec, EHc, EEc, C0c, C1c, C2c, C3c, F1c, F2c, F3c, SR2c, FRc;

static void tgskrg(double f, double d, double *pp, double *qq)
{   double cf, sf, ccf, t, p, q, g1, g2, g3, g4, g5, dd;
    cf = cos(f); sf = sin(f); ccf = cf*cf; t = (sf/cf)*(sf/cf);
    p = EEc*ccf/(1.-EEc); q = Ac / sqrt(1. - EEc*sf*sf);
    g1 = q*cf; g2 = g1*sf/2.; g3 = g1*ccf*(1.-t+p)/6.;
    g4 = g2*ccf*(5.-t)/12.; g5 = g1*ccf*ccf*(5.- 18.*t + t*t) / 120.;
    p = ((C2c - C3c*ccf)*ccf - C1c)*sf*cf + C0c*f;
    dd = d; q = g1*dd;
    dd *= d; p += g2*dd; dd *= d; q += g3*dd;
    dd *= d; p += g4*dd; dd *= d; q += g5*dd;
    *pp = p; *qq = q;
}

static void rgskrg(double p, double q, double *ff, double *dd)
{   double f, cf, sf, ccf, t, u, v, ec, tt, qq;
    f = p/C0c; cf = cos(f); sf = sin(f); ccf = cf*cf;
    f = ((F3c*ccf + F2c)*ccf + F1c)*sf*cf + f;
    cf = cos(f); sf = sin(f); t = sf/cf; u = 1. - EEc*sf*sf;
    v = t * u/(1.-EEc); ec = EEc*cf*cf/(1.-EEc); tt = t*t;
    q /= Ac/sqrt(u); qq = q*q;
    *ff = f - v*qq*(0.5 - qq*(5.+tt*3.)/24.);
    *dd = q*(1. - qq*((1.+tt+tt+ec) - qq*(5.+tt*28.+tt*tt*24.)/20.)/6.)/cf;
}

static void scvjtm(double gi, double gk)
{   double f, cf, sf, ccf;
    if (kw == 1) {
        Ac = AW; EEc = EEW; C0c = C0W; C1c = C1W; C2c = C2W; C3c = C3W;
        F1c = F1W; F2c = F2W; F3c = F3W;
    }
    else {
        Ac = AB; EEc = EEB; C0c = C0B; C1c = C1B; C2c = C2B; C3c = C3B;
        F1c = F1B; F2c = F2B; F3c = F3B;
    }
    f = gi/RAD; cf = cos(f); sf = sin(f); ccf = cf*cf;
    gx = ((C2c - C3c*ccf)*ccf - C1c)*sf*cf + C0c*f; glon = gk;
}

static void xcvjtm(double fi, double fk, double *x, double *y)
{   double p, q; tgskrg(fi/RAD, (fk-glon)/RAD, &p, &q);
    *x = (p-gx)*0.9999; *y = q*0.9999; }

static void zcvjtm(double x, double y, double *fi, double *fk)
{   double flat, dlon; rgskrg(x/0.9999+gx, y/0.9999, &flat, &dlon);
    *fi = flat*RAD; *fk = dlon*RAD+glon; }

static void scvutm(double gk)
{   if (kw == 1) {
        Ac = AW; EEc = EEW; C0c = C0W; C1c = C1W; C2c = C2W; C3c = C3W;
        F1c = F1W; F2c = F2W; F3c = F3W;
    }
    else {
        Ac = AB; EEc = EEB; C0c = C0B; C1c = C1B; C2c = C2B; C3c = C3B;
        F1c = F1B; F2c = F2B; F3c = F3B;
    }
    glon = gk;
}

static void xcvutm(double fi, double fk, double *x, double *y)
{   double p, q; tgskrg(fi/RAD, (fk-glon)/RAD, &p, &q);
    *x = p*0.9996; *y = q*0.9996 + 500.; }

static void zcvutm(double x, double y, double *fi, double *fk)
{   double flat, dlon; rgskrg(x/0.9996, (y-500.)/0.9996, &flat, &dlon);
    *fi = flat*RAD; *fk = dlon*RAD+glon; }

static void scvups(int ns)
{   double f, cf, sf, ccf;
    if (kw == 1) { Ac = AW; Ec = EW; EHc = EHW; }
    else { Ac = AB; Ec = EB; EHc = EHB; }
    ct = Ac*1.988/sqrt(pow((1.+Ec),(1.+Ec))*pow((1.-Ec),(1.-Ec)));
    st = (double)(3 - ns*2);
}

static void xcvups(double fi, double fk, double *x, double *y)
{   double cf, sf, f, es, t;
    cf = cos(fk/RAD); sf = sin(fk/RAD); f = fi/RAD;
    if (st > 0.) cf = -cf; else f = -f;
    es = Ec*sin(f); t = ct * tan((HPI-f)/2.) * pow(((1.+es)/(1.-es)),EHc);
    *x = t*cf + 2000.; *y = t*sf + 2000.;
}

static void zcvups(double x, double y, double *fi, double *fk)
{   double t0, t1, t2, es;
    x = x-2000.; y = y-2000.; t0 = sqrt(x*x+y*y)/ct;
    t2 = HPI - atan(t0)*2.;
    do{ t1 = t2; es = Ec*sin(t1);
        t2 = HPI - atan(t0*pow(((1.-es)/(1.+es)),EHc))*2.;
    } while (abs(t2-t1) > CR);
    if (st > 0.) { *fi = t2*RAD; *fk = atan2(y,-x)*RAD; }
    else         { *fi = -t2*RAD; *fk = atan2(y,x)*RAD; }
}

static void smercr(double gi, double gk)
{   double g, cg, sg, es, xp;
    if (kw == 1) { Ac = AW; Ec = EW; EHc = EHW; }
    else { Ac = AB; Ec = EB; EHc = EHB; }
    g = gi/RAD; cg = cos(g); es = Ec*sin(g); xp = (HPI+g)/2.;
    st = Ac*cg / sqrt(1.-es*es);
    gx = log(tan(xp) * pow(((1.-es)/(1.+es)),EHc)); glon = gk;
}

static void xmercr(double fi, double fk, double *x, double *y)
{   double f, es, xp, p;
    f = fi/RAD; es = Ec*sin(f); xp = (HPI+f)/2.;
    p = log(tan(xp) * pow(((1.-es)/(1.+es)),EHc));
    *x = (p-gx) * st; *y = (fk-glon)/RAD * st;
}

static void zmercr(double x, double y, double *fi, double *fk)
{   double t0, t1, t2, es;
    t0 = exp(x/st+gx); t2 = atan(t0)*2. - HPI;
    do{ t1 = t2; es = Ec*sin(t1);
        t2 = atan(t0*pow(((1.+es)/(1.-es)),EHc))*2. - HPI;
    } while (abs(t2-t1) > CR);
    *fi = t2*RAD; *fk = (y/st)*RAD + glon;
}

static void slccon(double gi, double gk, double sp1, double sp2)
{   double p0, es0, c0, t0, p1, es1, c1, t1, p2, es2, c2, t2;
    if (kw == 1) { Ac = AW; Ec = EW; EHc = EHW; }
    else { Ac = AB; Ec = EB; EHc = EHB; }
    p0 = gi/RAD; es0 = Ec*sin(p0);
    t0 = tan((HPI-p0)/2.)*pow(((1.+es0)/(1.-es0)),EHc);
    p1 = sp1/RAD; es1 = Ec*sin(p1); c1 = cos(p1)/sqrt(1.-es1*es1);
    t1 = tan((HPI-p1)/2.)*pow(((1.+es1)/(1.-es1)),EHc);
    if (sp1 == sp2) st = sin(sp1/RAD);
    else {
        p2 = sp2/RAD; es2 = Ec*sin(p2); c2 = cos(p2)/sqrt(1.-es2*es2);
        t2 = tan((HPI-p2)/2.)*pow(((1.+es2)/(1.-es2)),EHc);
        st = log(c1/c2)/log(t1/t2);
    }
    ct = (Ac/st)*c1/pow(t1,st);
    gx = ct*pow(t0,st); glon = gk;
}

static void xlccon(double fi, double fk, double *x, double *y)
{   double p, es, t, h, r;
    p = fi/RAD; es = Ec*sin(p); t = tan((HPI-p)/2.)*pow(((1.+es)/(1.-es)),EHc);
    h = st*(fk-glon)/RAD; r = ct*pow(t,st);
    *x = gx - r*cos(h); *y = r*sin(h);
}

static void zlccon(double x, double y, double *fi, double *fk)
{   double t0, t1, t2, es;
    x = gx-x; t0 = pow(sqrt(x*x+y*y)/ct,1./st);
    t2 = HPI - atan(t0)*2.;
    do{ t1 = t2; es = Ec*sin(t1);
        t2 = HPI - atan(t0*pow(((1.-es)/(1.+es)),EHc))*2.;
    } while (abs(t2-t1) > CR);
    *fi = t2*RAD; *fk = atan2(y,x)/st*RAD + glon;
}

static void slazea(double gi, double gk)
{   if (kw == 1) { SR2c = SR2W; FRc = FRW; } else { SR2c = SR2B; FRc = FRB; }
    ct = sin(gi/RAD); st = cos(gi/RAD); glon = gk; }

static void xlazea(double fi, double fk, double *x, double *y)
{   double cth, sth, dk, clon, slon, t, s;
    cth = sin(fi/RAD); sth = cos(fi/RAD);
    dk = fk-glon; clon = cos(dk/RAD); slon = sin(dk/RAD);
    t = 1. + ct*cth + st*sth*clon; s = SR2c / sqrt(t);
    *x = s * (st*cth - ct*sth*clon); *y = s * sth*slon;
}

static void zlazea(double x, double y, double *fi, double *fk)
{   double ss, s, cz, sz, cth, w;
    ss = x*x + y*y; s = sqrt(ss);
    cz = 1. - ss/(SR2c*SR2c); sz = sqrt(1.-cz*cz);
    if (ss == 0.) cth = ct; else cth = ct*cz + st*sz*x/s;
    ss = st*cz*s - ct*sz*x;
    *fi = asin(cth)*RAD; *fk = atan2(sz*y,ss)*RAD + glon;
}

/*----------------------------------------------------------------*/

int cvinit(int nca, double gi, double gk, double sp1, double sp2)
{   int n, m;
    if (nca < 200) { m = 0; n = nca; } else { m = 1; n = nca-200; }
    if ((n>0) && (n<=60)) { nc = n; kw = m; scvutm((nc*6-183)*60.); }
    else if (n ==  61) { nc = n; kw = m; scvups(1); }
    else if (n ==  62) { nc = n; kw = m; scvups(2); }
    else if (n ==  65) { nc = n; kw = m; scvutm(gk); }
    else if (n ==   0) { nc = n; kw = m; scvjtm(gi,gk); }
    else if (n ==  70) { nc = n; kw = m; smercr(gi,gk); }
    else if (n ==  71) { nc = n; kw = m; slccon(gi,gk,gi,gi); }
    else if (n ==  72) { nc = n; kw = m; slccon(gi,gk,sp1,sp2); }
    else if (n == 100) { nc = n; kw = m; slazea(gi,gk); }
    else if (n == 109) { nc = n; kw = m; slazea(gi,gk); }
    else if (n == 199) { nc = n; kw = m; gx = gi; glon = gk; }
    else { fprintf(stderr, "proj. #%d not implemented\n", nca); return(-1); }
    return(0);
}

int cviken(double fi, double fk, double *ye, double *xn)
{   int rc;
    if ((nc>0) &&(nc<=60)) xcvutm(fi,fk,xn,ye);
    else if (nc ==  61) xcvups(fi,fk,xn,ye);
    else if (nc ==  62) xcvups(fi,fk,xn,ye);
    else if (nc ==  65) xcvutm(fi,fk,xn,ye);
    else if (nc ==   0) xcvjtm(fi,fk,xn,ye);
    else if (nc ==  70) xmercr(fi,fk,xn,ye);
    else if (nc ==  71) xlccon(fi,fk,xn,ye);
    else if (nc ==  72) xlccon(fi,fk,xn,ye);
    else if (nc == 100) xlazea(fi,fk,xn,ye);
    else if (nc == 109) { xlazea(fi,fk,xn,ye); *ye /= FRc; *xn /= FRc; }
    else if (nc == 199) { *xn = fi-gx; *ye = fk-glon; }
    else { fprintf(stderr, "proj. coordinate undefined\n"); return(-1); }
    return(0);
}

int cvenik(double ye, double xn, double *fi, double *fk)
{   double x, y; int rc;
    if ((nc>0) && (nc<=60)) zcvutm(xn,ye,fi,fk);
    else if (nc ==  61) zcvups(xn,ye,fi,fk);
    else if (nc ==  62) zcvups(xn,ye,fi,fk);
    else if (nc ==  65) zcvutm(xn,ye,fi,fk);
    else if (nc ==   0) zcvjtm(xn,ye,fi,fk);
    else if (nc ==  70) zmercr(xn,ye,fi,fk);
    else if (nc ==  71) zlccon(xn,ye,fi,fk);
    else if (nc ==  72) zlccon(xn,ye,fi,fk);
    else if (nc == 100) zlazea(xn,ye,fi,fk);
    else if (nc == 109) zlazea(xn*FRc,ye*FRc,fi,fk);
    else if (nc == 199) { *fi = xn+gx; *fk = ye+glon; }
    else { fprintf(stderr, "proj. coordinate undefined\n"); return(-1); }
    return(0);
}

/*******************************************/
/*   FORTRAN Language Interface Routines   */
/*******************************************/

void cvinit_(int *nca, float *gi, float *gk, float *sp1, float *sp2)
{   int rc;
    if (((*nca>0) && (*nca<=62)) || ((*nca>200) && (*nca<=262)))
        rc = cvinit(*nca,0.,0.,0.,0.);
    else if ((*nca==72) || (*nca==272))
        rc = cvinit(*nca,(double)*gi,(double)*gk,(double)*sp1,(double)*sp2);
    else rc = cvinit(*nca,(double)*gi,(double)*gk,0.,0.);
    if (rc != 0) exit(1);
}

void cvdinit_(int *nca, double *gi, double *gk, double *sp1, double *sp2)
{   int rc;
    if (((*nca>0) && (*nca<=62)) || ((*nca>200) && (*nca<=262)))
        rc = cvinit(*nca,0.,0.,0.,0.);
    else if ((*nca==72) || (*nca==272)) rc = cvinit(*nca,*gi,*gk,*sp1,*sp2);
    else rc = cvinit(*nca,*gi,*gk,0.,0.);
    if (rc != 0) exit(1);
}

void cviken_(float *fi, float *fk, float *ye, float *xn)
{   double dye, dxn;
    if (cviken((double)*fi,(double)*fk,&dye,&dxn) != 0) exit(1);
    *ye = dye; *xn = dxn;
}

void cvdiken_(double *fi, double *fk, double *ye, double *xn)
{   if (cviken(*fi,*fk,ye,xn) != 0) exit(1); }

void cvenik_(float *ye, float *xn, float *fi, float *fk)
{   double di, dk;
    if (cvenik((double)*ye,(double)*xn,&di,&dk) != 0) exit(1);
    *fi = di; *fk = dk;
}

void cvdenik_(double *ye, double *xn, double *fi, double *fk)
{   if (cvenik(*ye,*xn,fi,fk) != 0) exit(1); }

/*----------------------------------------------------------------*/

void xyconv_(float *gi, float *gk, float *fi, float *fk, float *x, float *y)
{   double dx, dy;
    cvinit(0,(double)*gi,(double)*gk,0.,0.);
    cviken((double)*fi,(double)*fk,&dy,&dx);
    *x = dx; *y = dy;
}

void nxycnv_(float *gi, float *gk, float *fi, float *fk, float *x, float *y)
{   double dx, dy;
    cvinit(200,(double)*gi,(double)*gk,0.,0.);
    cviken((double)*fi,(double)*fk,&dy,&dx);
    *x = dx; *y = dy;
}

void utm_(int *nz, float *fi, float *fk, float *x, float *y)
{   double dx, dy;
    int n;
    if (*nz < 200) n = *nz; else n = *nz - 200;
    if ((n<1) || (n>60))
        { fprintf(stderr, "illegal zone (utm) : %d\n", n); exit(1); }
    cvinit(*nz,0.,0.,0.,0.);
    cviken((double)*fi,(double)*fk,&dy,&dx);
    *x = dx; *y = dy;
}

void ikconv_(float *gi, float *gk, float *x, float *y, float *fi, float *fk)
{   double di, dk;
    cvinit(0,(double)*gi,(double)*gk,0.,0.);
    cvenik((double)*y,(double)*x,&di,&dk);
    *fi = di; *fk = dk;
}

void nikcnv_(float *gi, float *gk, float *x, float *y, float *fi, float *fk)
{   double di, dk;
    cvinit(200,(double)*gi,(double)*gk,0.,0.);
    cvenik((double)*y,(double)*x,&di,&dk);
    *fi = di; *fk = dk;
}

void utmik_(int *nz, float *x, float *y, float *fi, float *fk)
{   double di, dk;
    int n;
    if (*nz < 200) n = *nz; else n = *nz - 200;
    if ((n<1) || (n>60))
        { fprintf(stderr, "illegal zone (utmik) : %d\n", n); exit(1); }
    cvinit(*nz,0.,0.,0.,0.);
    cvenik((double)*y,(double)*x,&di,&dk);
    *fi = di; *fk = dk;
}