/Proj4/PJ_krovak.c
C | 280 lines | 114 code | 57 blank | 109 comment | 8 complexity | b535e1b3d33c62a0045281ee4e163810 MD5 | raw file
1/****************************************************************************** 2 * $Id: PJ_krovak.c,v 1.9 2007/03/07 17:32:32 fwarmerdam Exp $ 3 * 4 * Project: PROJ.4 5 * Purpose: Implementation of the krovak (Krovak) projection. 6 * Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3 7 * Author: Thomas Flemming, tf@ttqv.com 8 * 9 ****************************************************************************** 10 * Copyright (c) 2001, Thomas Flemming, tf@ttqv.com 11 * 12 * Permission is hereby granted, free of charge, to any person obtaining a 13 * copy of this software and associated documentation files (the "Software"), 14 * to deal in the Software without restriction, including without limitation 15 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 16 * and/or sell copies of the Software, and to permit persons to whom the 17 * Software is furnished to do so, subject to the following conditions: 18 * 19 * The above copyright notice and this permission notice shall be included 20 * in all copies or substantial portions of the Software. 21 * 22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 23 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 24 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 25 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS 26 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN 27 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN 28 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 29 * SOFTWARE. 30 ****************************************************************************** 31 * 32 * $Log: PJ_krovak.c,v $ 33 * Revision 1.9 2007/03/07 17:32:32 fwarmerdam 34 * remove orphan semicolon. 35 * 36 * Revision 1.8 2007/03/07 17:28:08 fwarmerdam 37 * Make it reasonably clear that this is ellipsoidal in the code. 38 * 39 * Revision 1.7 2007/03/07 17:25:34 fwarmerdam 40 * report krovak as ellipsoidal, not spherical 41 * 42 * Revision 1.6 2006/09/14 13:10:50 fwarmerdam 43 * Add +czech flag to control reversal of signs (bug 1133,147) 44 * 45 * Revision 1.5 2006/03/30 01:22:48 fwarmerdam 46 * Removed win32 only debug hack. 47 * 48 * Revision 1.4 2002/12/15 22:31:04 warmerda 49 * handle lon_0, k, and prime meridian properly 50 * 51 * Revision 1.3 2002/12/15 00:13:30 warmerda 52 * lat_0 may now be set by user, but still defaults to 49d30N 53 * 54 * Revision 1.2 2002/12/14 19:35:21 warmerda 55 * updated headers 56 * 57 */ 58 59#define PROJ_PARMS__ \ 60 double C_x; 61#define PJ_LIB__ 62 63#include "projects.h" 64#include <string.h> 65#include <stdio.h> 66 67PJ_CVSID("$Id: PJ_krovak.c,v 1.9 2007/03/07 17:32:32 fwarmerdam Exp $"); 68 69PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Ellps."; 70 71/** 72 NOTES: According to EPSG the full Krovak projection method should have 73 the following parameters. Within PROJ.4 the azimuth, and pseudo 74 standard parallel are hardcoded in the algorithm and can't be 75 altered from outside. The others all have defaults to match the 76 common usage with Krovak projection. 77 78 lat_0 = latitude of centre of the projection 79 80 lon_0 = longitude of centre of the projection 81 82 ** = azimuth (true) of the centre line passing through the centre of the projection 83 84 ** = latitude of pseudo standard parallel 85 86 k = scale factor on the pseudo standard parallel 87 88 x_0 = False Easting of the centre of the projection at the apex of the cone 89 90 y_0 = False Northing of the centre of the projection at the apex of the cone 91 92 **/ 93 94 95 96FORWARD(e_forward); /* ellipsoid */ 97/* calculate xy from lat/lon */ 98 99 char errmess[255]; 100 char tmp[16]; 101 102/* Constants, identical to inverse transform function */ 103 double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n; 104 double gfi, u, fi0, lon17, lamdd, deltav, s, d, eps, ro; 105 106 107 s45 = 0.785398163397448; /* 45� */ 108 s90 = 2 * s45; 109 fi0 = P->phi0; /* Latitude of projection centre 49� 30' */ 110 111 /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 112 be set to 1 here. 113 Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128, 114 e2=0.006674372230614; 115 */ 116 a = 1; /* 6377397.155; */ 117 /* e2 = P->es;*/ /* 0.006674372230614; */ 118 e2 = 0.006674372230614; 119 e = sqrt(e2); 120 121 alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2)); 122 123 uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ 124 u0 = asin(sin(fi0) / alfa); 125 g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. ); 126 127 k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g; 128 129 k1 = P->k0; 130 n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2)); 131 s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78� 30'00" N */ 132 n = sin(s0); 133 ro0 = k1 * n0 / tan(s0); 134 ad = s90 - uq; 135 136/* Transformation */ 137 138 gfi =pow ( ((1. + e * sin(lp.phi)) / 139 (1. - e * sin(lp.phi))) , (alfa * e / 2.)); 140 141 u= 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45); 142 143 deltav = - lp.lam * alfa; 144 145 s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav)); 146 d = asin(cos(u) * sin(deltav) / cos(s)); 147 eps = n * d; 148 ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n) ; 149 150 /* x and y are reverted! */ 151 xy.y = ro * cos(eps) / a; 152 xy.x = ro * sin(eps) / a; 153 154 if( !pj_param(P -> params, "tczech").i ) 155 { 156 xy.y *= -1.0; 157 xy.x *= -1.0; 158 } 159 160 return (xy); 161} 162 163 164 165INVERSE(e_inverse); /* ellipsoid */ 166 /* calculate lat/lon from xy */ 167 168/* Constants, identisch wie in der Umkehrfunktion */ 169 double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n; 170 double u, l24, lamdd, deltav, s, d, eps, ro, fi1, xy0, lon17; 171 int ok; 172 173 s45 = 0.785398163397448; /* 45� */ 174 s90 = 2 * s45; 175 fi0 = P->phi0; /* Latitude of projection centre 49� 30' */ 176 177 178 /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 179 be set to 1 here. 180 Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128, 181 e2=0.006674372230614; 182 */ 183 a = 1; /* 6377397.155; */ 184 /* e2 = P->es; */ /* 0.006674372230614; */ 185 e2 = 0.006674372230614; 186 e = sqrt(e2); 187 188 alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2)); 189 uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ 190 u0 = asin(sin(fi0) / alfa); 191 g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. ); 192 193 k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g; 194 195 k1 = P->k0; 196 n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2)); 197 s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78� 30'00" N */ 198 n = sin(s0); 199 ro0 = k1 * n0 / tan(s0); 200 ad = s90 - uq; 201 202 203/* Transformation */ 204 /* revert y, x*/ 205 xy0=xy.x; 206 xy.x=xy.y; 207 xy.y=xy0; 208 209 if( !pj_param(P -> params, "tczech").i ) 210 { 211 xy.x *= -1.0; 212 xy.y *= -1.0; 213 } 214 215 ro = sqrt(xy.x * xy.x + xy.y * xy.y); 216 eps = atan2(xy.y, xy.x); 217 d = eps / sin(s0); 218 s = 2. * (atan( pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45); 219 220 u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d)); 221 deltav = asin(cos(s) * sin(d) / cos(u)); 222 223 lp.lam = P->lam0 - deltav / alfa; 224 225/* ITERATION FOR lp.phi */ 226 fi1 = u; 227 228 ok = 0; 229 do 230 { 231 lp.phi = 2. * ( atan( pow( k, -1. / alfa) * 232 pow( tan(u / 2. + s45) , 1. / alfa) * 233 pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.) 234 ) - s45); 235 236 if (fabs(fi1 - lp.phi) < 0.000000000000001) ok=1; 237 fi1 = lp.phi; 238 239 } 240 while (ok==0); 241 242 lp.lam -= P->lam0; 243 244 return (lp); 245} 246 247FREEUP; if (P) pj_dalloc(P); } 248 249ENTRY0(krovak) 250 double ts; 251 /* read some Parameters, 252 * here Latitude Truescale */ 253 254 ts = pj_param(P->params, "rlat_ts").f; 255 P->C_x = ts; 256 257 /* we want Bessel as fixed ellipsoid */ 258 P->a = 6377397.155; 259 P->e = sqrt(P->es = 0.006674372230614); 260 261 /* if latitude of projection center is not set, use 49d30'N */ 262 if (!pj_param(P->params, "tlat_0").i) 263 P->phi0 = 0.863937979737193; 264 265 /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */ 266 /* that will correspond to using longitudes relative to greenwich */ 267 /* as input and output, instead of lat/long relative to Ferro */ 268 if (!pj_param(P->params, "tlon_0").i) 269 P->lam0 = 0.7417649320975901 - 0.308341501185665; 270 271 /* if scale not set default to 0.9999 */ 272 if (!pj_param(P->params, "tk").i) 273 P->k0 = 0.9999; 274 275 /* always the same */ 276 P->inv = e_inverse; 277 P->fwd = e_forward; 278 279ENDENTRY(P) 280