~ubuntu-branches/ubuntu/vivid/proj/vivid

« back to all changes in this revision

Viewing changes to src/PJ_lsat.c

Committer: Bazaar Package Importer
Author(s): Peter S Galbraith
Date: 2002-01-11 10:27:12 UTC
Revision ID: james.westby@ubuntu.com-20020111102712-ayi18r8y2eesv0y9

Tags: upstream-4.4.5

Import upstream version 4.4.5

files added:

AUTHORS

COPYING

ChangeLog

INSTALL

Makefile.am

Makefile.in

NEWS

README

aclocal.m4

config.guess

config.sub

configure

configure.in

install-sh

ltconfig

ltmain.sh

man/Makefile.am

man/Makefile.in

man/man1

man/man1/Makefile.am

man/man1/Makefile.in

man/man1/cs2cs.1

man/man1/geod.1

man/man1/nad2nad.1

man/man1/proj.1

man/man3

man/man3/Makefile.am

man/man3/Makefile.in

man/man3/pj_init.3

missing

mkinstalldirs

nad/GL27

nad/Makefile.am

nad/Makefile.in

nad/README

nad/epsg

nad/nad.lst

nad/nad27

nad/nad83

nad/pj_out27.dist

nad/pj_out83.dist

nad/proj_def.dat

nad/test27

nad/test83

nad/world

src/Makefile.am

src/Makefile.in

src/PJ_aea.c

src/PJ_aeqd.c

src/PJ_airy.c

src/PJ_aitoff.c

src/PJ_august.c

src/PJ_bacon.c

src/PJ_bipc.c

src/PJ_boggs.c

src/PJ_bonne.c

src/PJ_cass.c

src/PJ_cc.c

src/PJ_cea.c

src/PJ_chamb.c

src/PJ_collg.c

src/PJ_crast.c

src/PJ_denoy.c

src/PJ_eck1.c

src/PJ_eck2.c

src/PJ_eck3.c

src/PJ_eck4.c

src/PJ_eck5.c

src/PJ_eqc.c

src/PJ_eqdc.c

src/PJ_fahey.c

src/PJ_fouc_s.c

src/PJ_gall.c

src/PJ_gins8.c

src/PJ_gn_sinu.c

src/PJ_gnom.c

src/PJ_goode.c

src/PJ_hammer.c

src/PJ_hatano.c

src/PJ_imw_p.c

src/PJ_labrd.c

src/PJ_laea.c

src/PJ_lagrng.c

src/PJ_larr.c

src/PJ_lask.c

src/PJ_lcc.c

src/PJ_loxim.c

src/PJ_lsat.c

src/PJ_mbt_fps.c

src/PJ_mbtfpp.c

src/PJ_mbtfpq.c

src/PJ_merc.c

src/PJ_mill.c

src/PJ_mod_ster.c

src/PJ_moll.c

src/PJ_mpoly.c

src/PJ_nell.c

src/PJ_nell_h.c

src/PJ_nocol.c

src/PJ_nsper.c

src/PJ_nzmg.c

src/PJ_ob_tran.c

src/PJ_ocea.c

src/PJ_oea.c

src/PJ_omerc.c

src/PJ_ortho.c

src/PJ_poly.c

src/PJ_putp2.c

src/PJ_putp3.c

src/PJ_putp4p.c

src/PJ_putp5.c

src/PJ_putp6.c

src/PJ_robin.c

src/PJ_rpoly.c

src/PJ_sconics.c

src/PJ_somerc.c

src/PJ_stere.c

src/PJ_sts.c

src/PJ_tcc.c

src/PJ_tcea.c

src/PJ_tmerc.c

src/PJ_tpeqd.c

src/PJ_urm5.c

src/PJ_urmfps.c

src/PJ_vandg.c

src/PJ_vandg2.c

src/PJ_vandg4.c

src/PJ_wag2.c

src/PJ_wag3.c

src/PJ_wag7.c

src/PJ_wink1.c

src/PJ_wink2.c

src/aasincos.c

src/adjlon.c

src/bch2bps.c

src/bchgen.c

src/biveval.c

src/cs2cs.c

src/dmstor.c

src/emess.c

src/emess.h

src/gen_cheb.c

src/geocent.c

src/geocent.h

src/geod.c

src/geod_for.c

src/geod_inv.c

src/geod_set.c

src/geodesic.h

src/makefile.vc

src/mk_cheby.c

src/nad2bin.c

src/nad2nad.c

src/nad_cvt.c

src/nad_init.c

src/nad_intr.c

src/nad_list.h

src/p_series.c

src/pj_apply_gridshift.c

src/pj_auth.c

src/pj_datum_set.c

src/pj_datums.c

src/pj_deriv.c

src/pj_ell_set.c

src/pj_ellps.c

src/pj_errno.c

src/pj_factors.c

src/pj_fwd.c

src/pj_init.c

src/pj_inv.c

src/pj_latlong.c

src/pj_list.c

src/pj_list.h

src/pj_malloc.c

src/pj_mlfn.c

src/pj_msfn.c

src/pj_open_lib.c

src/pj_param.c

src/pj_phi2.c

src/pj_pr_list.c

src/pj_qsfn.c

src/pj_release.c

src/pj_strerrno.c

src/pj_transform.c

src/pj_tsfn.c

src/pj_units.c

src/pj_utils.c

src/pj_zpoly1.c

src/proj.c

src/proj.def

src/proj_api.h

src/proj_config.h.in

src/projects.h

src/rtodms.c

src/stamp-h.in

src/strtod.c

src/vector1.c

Show diffs side-by-side

added added

removed removed

src/PJ_lsat.c

#ifndef lint

static const char SCCSID[]="@(#)PJ_lsat.c 4.1 94/02/15 GIE REL";

#endif

/* based upon Snyder and Linck, USGS-NMD */

#define PROJ_PARMS__ \

double a2, a4, b, c1, c3; \

double q, t, u, w, p22, sa, ca, xj, rlm, rlm2;

#define PJ_LIB__

#include <projects.h>

PROJ_HEAD(lsat, "Space oblique for LANDSAT")

"\n\tCyl, Sph&Ell\n\tlsat= path=";

#define TOL 1e-7

#define PI_HALFPI 4.71238898038468985766

#define TWOPI_HALFPI 7.85398163397448309610

static void

seraz0(double lam, double mult, PJ *P) {

double sdsq, h, s, fc, sd, sq, d__1;

lam *= DEG_TO_RAD;

sd = sin(lam);

sdsq = sd * sd;

s = P->p22 * P->sa * cos(lam) * sqrt((1. + P->t * sdsq) / ((

1. + P->w * sdsq) * (1. + P->q * sdsq)));

d__1 = 1. + P->q * sdsq;

h = sqrt((1. + P->q * sdsq) / (1. + P->w * sdsq)) * ((1. +

P->w * sdsq) / (d__1 * d__1) - P->p22 * P->ca);

sq = sqrt(P->xj * P->xj + s * s);

P->b += fc = mult * (h * P->xj - s * s) / sq;

P->a2 += fc * cos(lam + lam);

P->a4 += fc * cos(lam * 4.);

fc = mult * s * (h + P->xj) / sq;

P->c1 += fc * cos(lam);

P->c3 += fc * cos(lam * 3.);

}

FORWARD(e_forward); /* ellipsoid */

int l, nn;

double lamt, xlam, sdsq, c, d, s, lamdp, phidp, lampp, tanph,

lamtp, cl, sd, sp, fac, sav, tanphi;

if (lp.phi > HALFPI)

lp.phi = HALFPI;

else if (lp.phi < -HALFPI)

lp.phi = -HALFPI;

lampp = lp.phi >= 0. ? HALFPI : PI_HALFPI;

tanphi = tan(lp.phi);

for (nn = 0;;) {

sav = lampp;

lamtp = lp.lam + P->p22 * lampp;

cl = cos(lamtp);

if (fabs(cl) < TOL)

lamtp -= TOL;

fac = lampp - sin(lampp) * (cl < 0. ? -HALFPI : HALFPI);

for (l = 50; l; --l) {

lamt = lp.lam + P->p22 * sav;

if (fabs(c = cos(lamt)) < TOL)

lamt -= TOL;

xlam = (P->one_es * tanphi * P->sa + sin(lamt) * P->ca) / c;

lamdp = atan(xlam) + fac;

if (fabs(fabs(sav) - fabs(lamdp)) < TOL)

break;

sav = lamdp;

}

if (!l || ++nn >= 3 || (lamdp > P->rlm && lamdp < P->rlm2))

break;

if (lamdp <= P->rlm)

lampp = TWOPI_HALFPI;

else if (lamdp >= P->rlm2)

lampp = HALFPI;

}

if (l) {

sp = sin(lp.phi);

phidp = aasin((P->one_es * P->ca * sp - P->sa * cos(lp.phi) *

sin(lamt)) / sqrt(1. - P->es * sp * sp));

tanph = log(tan(FORTPI + .5 * phidp));

sd = sin(lamdp);

sdsq = sd * sd;

s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)

/ ((1. + P->w * sdsq) * (1. + P->q * sdsq)));

d = sqrt(P->xj * P->xj + s * s);

xy.x = P->b * lamdp + P->a2 * sin(2. * lamdp) + P->a4 *

sin(lamdp * 4.) - tanph * s / d;

xy.y = P->c1 * sd + P->c3 * sin(lamdp * 3.) + tanph * P->xj / d;

} else

xy.x = xy.y = HUGE_VAL;

return xy;

}

INVERSE(e_inverse); /* ellipsoid */

int nn;

double lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;

lamdp = xy.x / P->b;

nn = 50;

do {

sav = lamdp;

sd = sin(lamdp);

sdsq = sd * sd;

s = P->p22 * P->sa * cos(lamdp) * sqrt((1. + P->t * sdsq)

/ ((1. + P->w * sdsq) * (1. + P->q * sdsq)));

lamdp = xy.x + xy.y * s / P->xj - P->a2 * sin(

100

2. * lamdp) - P->a4 * sin(lamdp * 4.) - s / P->xj * (

101

P->c1 * sin(lamdp) + P->c3 * sin(lamdp * 3.));

102

lamdp /= P->b;

103

} while (fabs(lamdp - sav) >= TOL && --nn);

104

sl = sin(lamdp);

105

fac = exp(sqrt(1. + s * s / P->xj / P->xj) * (xy.y -

106

P->c1 * sl - P->c3 * sin(lamdp * 3.)));

107

phidp = 2. * (atan(fac) - FORTPI);

108

dd = sl * sl;

109

if (fabs(cos(lamdp)) < TOL)

110

lamdp -= TOL;

111

spp = sin(phidp);

112

sppsq = spp * spp;

113

lamt = atan(((1. - sppsq * P->rone_es) * tan(lamdp) *

114

P->ca - spp * P->sa * sqrt((1. + P->q * dd) * (

115

1. - sppsq) - sppsq * P->u) / cos(lamdp)) / (1. - sppsq

116

* (1. + P->u)));

117

sl = lamt >= 0. ? 1. : -1.;

118

scl = cos(lamdp) >= 0. ? 1. : -1;

119

lamt -= HALFPI * (1. - scl) * sl;

120

lp.lam = lamt - P->p22 * lamdp;

121

if (fabs(P->sa) < TOL)

122

lp.phi = aasin(spp / sqrt(P->one_es * P->one_es + P->es * sppsq));

123

else

124

lp.phi = atan((tan(lamdp) * cos(lamt) - P->ca * sin(lamt)) /

125

(P->one_es * P->sa));

126

return lp;

127

}

128

FREEUP; if (P) pj_dalloc(P); }

129

ENTRY0(lsat)

130

int land, path;

131

double lam, alf, esc, ess;

132

133

land = pj_param(P->params, "ilsat").i;

134

if (land <= 0 || land > 5) E_ERROR(-28);

135

path = pj_param(P->params, "ipath").i;

136

if (path <= 0 || path > (land <= 3 ? 251 : 233)) E_ERROR(-29);

137

if (land <= 3) {

138

P->lam0 = DEG_TO_RAD * 128.87 - TWOPI / 251. * path;

139

P->p22 = 103.2669323;

140

alf = DEG_TO_RAD * 99.092;

141

} else {

142

P->lam0 = DEG_TO_RAD * 129.3 - TWOPI / 233. * path;

143

P->p22 = 98.8841202;

144

alf = DEG_TO_RAD * 98.2;

145

}

146

P->p22 /= 1440.;

147

P->sa = sin(alf);

148

P->ca = cos(alf);

149

if (fabs(P->ca) < 1e-9)

150

P->ca = 1e-9;

151

esc = P->es * P->ca * P->ca;

152

ess = P->es * P->sa * P->sa;

153

P->w = (1. - esc) * P->rone_es;

154

P->w = P->w * P->w - 1.;

155

P->q = ess * P->rone_es;

156

P->t = ess * (2. - P->es) * P->rone_es * P->rone_es;

157

P->u = esc * P->rone_es;

158

P->xj = P->one_es * P->one_es * P->one_es;

159

P->rlm = PI * (1. / 248. + .5161290322580645);

160

P->rlm2 = P->rlm + TWOPI;

161

P->a2 = P->a4 = P->b = P->c1 = P->c3 = 0.;

162

seraz0(0., 1., P);

163

for (lam = 9.; lam <= 81.0001; lam += 18.)

164

seraz0(lam, 4., P);

165

for (lam = 18; lam <= 72.0001; lam += 18.)

166

seraz0(lam, 2., P);

167

seraz0(90., 1., P);

168

P->a2 /= 30.;

169

P->a4 /= 60.;

170

P->b /= 30.;

171

P->c1 /= 15.;

172

P->c3 /= 45.;

173

P->inv = e_inverse; P->fwd = e_forward;

174

ENDENTRY(P)

Older »