Subversion Repositories eFlore/Projets.eflore-projets

Rev

Rev 1001 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1001 delphine 1
<?php
2
 
3
/**
4
 * Author : Julien Moquet
5
 *
6
 * Inspired by Proj4js from Mike Adair madairATdmsolutions.ca
7
 *                      and Richard Greenwood rich@greenwoodmap.com
8
 * License: LGPL as per: http://www.gnu.org/copyleft/lesser.html
9
 */
10
class Proj4phpCommon {
11
 
12
    public $PI = M_PI; #3.141592653589793238; //Math.PI,
13
    public $HALF_PI = M_PI_2; #1.570796326794896619; //Math.PI*0.5,
14
    public $TWO_PI = 6.283185307179586477; //Math.PI*2,
15
    public $FORTPI = 0.78539816339744833;
16
    public $R2D = 57.29577951308232088;
17
    public $D2R = 0.01745329251994329577;
18
    public $SEC_TO_RAD = 4.84813681109535993589914102357e-6; /* SEC_TO_RAD = Pi/180/3600 */
19
    public $EPSLN = 1.0e-10;
20
    public $MAX_ITER = 20;
21
    // following constants from geocent.c
22
    public $COS_67P5 = 0.38268343236508977;  /* cosine of 67.5 degrees */
23
    public $AD_C = 1.0026000;                /* Toms region 1 constant */
24
 
25
    /* datum_type values */
26
    public $PJD_UNKNOWN = 0;
27
    public $PJD_3PARAM = 1;
28
    public $PJD_7PARAM = 2;
29
    public $PJD_GRIDSHIFT = 3;
30
    public $PJD_WGS84 = 4;   // WGS84 or equivalent
31
    public $PJD_NODATUM = 5;   // WGS84 or equivalent
32
 
33
    const SRS_WGS84_SEMIMAJOR = 6378137.0;  // only used in grid shift transforms
34
 
35
    // ellipoid pj_set_ell.c
36
 
37
    public $SIXTH = .1666666666666666667; /* 1/6 */
38
    public $RA4 = .04722222222222222222; /* 17/360 */
39
    public $RA6 = .02215608465608465608; /* 67/3024 */
40
    public $RV4 = .06944444444444444444; /* 5/72 */
41
    public $RV6 = .04243827160493827160; /* 55/1296 */
42
 
43
 
44
    /* meridinal distance for ellipsoid and inverse
45
     * *	8th degree - accurate to < 1e-5 meters when used in conjuction
46
     * *		with typical major axis values.
47
     * *	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
48
     */
49
    protected $C00 = 1.0;
50
    protected $C02 = .25;
51
    protected $C04 = .046875;
52
    protected $C06 = .01953125;
53
    protected $C08 = .01068115234375;
54
    protected $C22 = .75;
55
    protected $C44 = .46875;
56
    protected $C46 = .01302083333333333333;
57
    protected $C48 = .00712076822916666666;
58
    protected $C66 = .36458333333333333333;
59
    protected $C68 = .00569661458333333333;
60
    protected $C88 = .3076171875;
61
 
62
    /**
63
     * Function to compute the constant small m which is the radius of
64
     * a parallel of latitude, phi, divided by the semimajor axis.
65
     *
66
     * @param type $eccent
67
     * @param type $sinphi
68
     * @param type $cosphi
69
     * @return type
70
     */
71
    public function msfnz( $eccent, $sinphi, $cosphi ) {
72
 
73
        $con = $eccent * $sinphi;
74
 
75
        return $cosphi / (sqrt( 1.0 - $con * $con ));
76
    }
77
 
78
    /**
79
     * Function to compute the constant small t for use in the forward
80
     * computations in the Lambert Conformal Conic and the Polar
81
     * Stereographic projections.
82
     *
83
     * @param type $eccent
84
     * @param type $phi
85
     * @param type $sinphi
86
     * @return type
87
     */
88
    public function tsfnz( $eccent, $phi, $sinphi ) {
89
 
90
        $con = $eccent * $sinphi;
91
        $com = 0.5 * $eccent;
92
        $con = pow( ((1.0 - $con) / (1.0 + $con) ), $com );
93
 
94
        return (tan( .5 * (M_PI_2 - $phi) ) / $con);
95
    }
96
 
97
    /**
98
     * Function to compute the latitude angle, phi2, for the inverse of the
99
     * Lambert Conformal Conic and Polar Stereographic projections.
100
     *
101
     * rise up an assertion if there is no convergence.
102
     *
103
     * @param type $eccent
104
     * @param type $ts
105
     * @return type
106
     */
107
    public function phi2z( $eccent, $ts ) {
108
 
109
        $eccnth = .5 * $eccent;
110
        $phi = M_PI_2 - 2 * atan( $ts );
111
 
112
        for( $i = 0; $i <= 15; $i++ ) {
113
            $con = $eccent * sin( $phi );
114
            $dphi = M_PI_2 - 2 * atan( $ts * (pow( ((1.0 - $con) / (1.0 + $con) ), $eccnth )) ) - $phi;
115
            $phi += $dphi;
116
            if( abs( $dphi ) <= .0000000001 )
117
                return $phi;
118
        }
119
        assert( "false; /* phi2z has NoConvergence */" );
120
 
121
        return (-9999);
122
    }
123
 
124
    /**
125
     * Function to compute constant small q which is the radius of a
126
     * parallel of latitude, phi, divided by the semimajor axis.
127
     *
128
     * @param type $eccent
129
     * @param type $sinphi
130
     * @return type
131
     */
132
    public function qsfnz( $eccent, $sinphi ) {
133
 
134
        if( $eccent > 1.0e-7 ) {
135
 
136
            $con = $eccent * $sinphi;
137
 
138
            return (( 1.0 - $eccent * $eccent) * ($sinphi / (1.0 - $con * $con) - (.5 / $eccent) * log( (1.0 - $con) / (1.0 + $con) )));
139
        }
140
 
141
        return (2.0 * $sinphi);
142
    }
143
 
144
    /**
145
     * Function to eliminate roundoff errors in asin
146
     *
147
     * @param type $x
148
     * @return type
149
     */
150
    public function asinz( $x ) {
151
 
152
        return asin(
153
            abs( $x ) > 1.0 ? ($x > 1.0 ? 1.0 : -1.0) : $x
154
        );
155
 
156
        #if( abs( $x ) > 1.0 ) {
157
        #    $x = ($x > 1.0) ? 1.0 : -1.0;
158
        #}
159
        #return asin( $x );
160
    }
161
 
162
    /**
163
     * following functions from gctpc cproj.c for transverse mercator projections
164
     *
165
     * @param type $x
166
     * @return type
167
     */
168
    public function e0fn( $x ) {
169
        return (1.0 - 0.25 * $x * (1.0 + $x / 16.0 * (3.0 + 1.25 * $x)));
170
    }
171
 
172
    /**
173
     *
174
     * @param type $x
175
     * @return type
176
     */
177
    public function e1fn( $x ) {
178
        return (0.375 * $x * (1.0 + 0.25 * $x * (1.0 + 0.46875 * $x)));
179
    }
180
 
181
    /**
182
     *
183
     * @param type $x
184
     * @return type
185
     */
186
    public function e2fn( $x ) {
187
        return (0.05859375 * $x * $x * (1.0 + 0.75 * $x));
188
    }
189
 
190
    /**
191
     *
192
     * @param type $x
193
     * @return type
194
     */
195
    public function e3fn( $x ) {
196
        return ($x * $x * $x * (35.0 / 3072.0));
197
    }
198
 
199
    /**
200
     *
201
     * @param type $e0
202
     * @param type $e1
203
     * @param type $e2
204
     * @param type $e3
205
     * @param type $phi
206
     * @return type
207
     */
208
    public function mlfn( $e0, $e1, $e2, $e3, $phi ) {
209
        return ($e0 * $phi - $e1 * sin( 2.0 * $phi ) + $e2 * sin( 4.0 * $phi ) - $e3 * sin( 6.0 * $phi ));
210
    }
211
 
212
    /**
213
     *
214
     * @param type $esinp
215
     * @param type $exp
216
     * @return type
217
     */
218
    public function srat( $esinp, $exp ) {
219
        return (pow( (1.0 - $esinp) / (1.0 + $esinp), $exp ));
220
    }
221
 
222
    /**
223
     * Function to return the sign of an argument
224
     *
225
     * @param type $x
226
     * @return type
227
     */
228
    public function sign( $x ) {
229
 
230
        return $x < 0.0 ? -1 : 1;
231
    }
232
 
233
    /**
234
     * Function to adjust longitude to -180 to 180; input in radians
235
     *
236
     * @param type $x
237
     * @return type
238
     */
239
    public function adjust_lon( $x ) {
240
 
241
        return (abs( $x ) < M_PI) ? $x : ($x - ($this->sign( $x ) * $this->TWO_PI) );
242
    }
243
 
244
    /**
245
     * IGNF - DGR : algorithms used by IGN France
246
     * Function to adjust latitude to -90 to 90; input in radians
247
     *
248
     * @param type $x
249
     * @return type
250
     */
251
    public function adjust_lat( $x ) {
252
 
253
        $x = (abs( $x ) < M_PI_2) ? $x : ($x - ($this->sign( $x ) * M_PI) );
254
 
255
        return $x;
256
    }
257
 
258
    /**
259
     * Latitude Isometrique - close to tsfnz ...
260
     *
261
     * @param type $eccent
262
     * @param float $phi
263
     * @param type $sinphi
264
     * @return string
265
     */
266
    public function latiso( $eccent, $phi, $sinphi ) {
267
 
268
        if( abs( $phi ) > M_PI_2 )
269
            return +NaN;
270
        if( $phi == M_PI_2 )
271
            return INF;
272
        if( $phi == -1.0 * M_PI_2 )
273
            return -1.0 * INF;
274
 
275
        $con = $eccent * $sinphi;
276
 
277
        return log( tan( (M_PI_2 + $phi) / 2.0 ) ) + $eccent * log( (1.0 - $con) / (1.0 + $con) ) / 2.0;
278
    }
279
 
280
    /**
281
     *
282
     * @param type $x
283
     * @param type $L
284
     * @return type
285
     */
286
    public function fL( $x, $L ) {
287
        return 2.0 * atan( $x * exp( $L ) ) - M_PI_2;
288
    }
289
 
290
    /**
291
     * Inverse Latitude Isometrique - close to ph2z
292
     *
293
     * @param type $eccent
294
     * @param type $ts
295
     * @return type
296
     */
297
    public function invlatiso( $eccent, $ts ) {
298
 
299
        $phi = $this->fL( 1.0, $ts );
300
        $Iphi = 0.0;
301
        $con = 0.0;
302
 
303
        do {
304
            $Iphi = $phi;
305
            $con = $eccent * sin( $Iphi );
306
            $phi = $this->fL( exp( $eccent * log( (1.0 + $con) / (1.0 - $con) ) / 2.0 ), $ts );
307
        } while( abs( $phi - $Iphi ) > 1.0e-12 );
308
 
309
        return $phi;
310
    }
311
 
312
    /**
313
     * Grande Normale
314
     *
315
     * @param type $a
316
     * @param type $e
317
     * @param type $sinphi
318
     * @return type
319
     */
320
    public function gN( $a, $e, $sinphi ) {
321
        $temp = $e * $sinphi;
322
        return $a / sqrt( 1.0 - $temp * $temp );
323
    }
324
 
325
    /**
326
     * code from the PROJ.4 pj_mlfn.c file;  this may be useful for other projections
327
     *
328
     * @param type $es
329
     * @return type
330
     */
331
    public function pj_enfn( $es ) {
332
 
333
        $en = array( );
334
        $en[0] = $this->C00 - $es * ($this->C02 + $es * ($this->C04 + $es * ($this->C06 + $es * $this->C08)));
335
        $en[1] = es * ($this->C22 - $es * ($this->C04 + $es * ($this->C06 + $es * $this->C08)));
336
        $t = $es * $es;
337
        $en[2] = $t * ($this->C44 - $es * ($this->C46 + $es * $this->C48));
338
        $t *= $es;
339
        $en[3] = $t * ($this->C66 - $es * $this->C68);
340
        $en[4] = $t * $es * $this->C88;
341
 
342
        return $en;
343
    }
344
 
345
    /**
346
     *
347
     * @param type $phi
348
     * @param type $sphi
349
     * @param type $cphi
350
     * @param type $en
351
     * @return type
352
     */
353
    public function pj_mlfn( $phi, $sphi, $cphi, $en ) {
354
 
355
        $cphi *= $sphi;
356
        $sphi *= $sphi;
357
 
358
        return ($en[0] * $phi - $cphi * ($en[1] + $sphi * ($en[2] + $sphi * ($en[3] + $sphi * $en[4]))));
359
    }
360
 
361
    /**
362
     *
363
     * @param type $arg
364
     * @param type $es
365
     * @param type $en
366
     * @return type
367
     */
368
    public function pj_inv_mlfn( $arg, $es, $en ) {
369
 
370
        $k = (float) 1 / (1 - $es);
371
        $phi = $arg;
372
 
373
        for( $i = Proj4php::$common->MAX_ITER; $i; --$i ) { /* rarely goes over 2 iterations */
374
            $s = sin( $phi );
375
            $t = 1. - $es * $s * $s;
376
            //$t = $this->pj_mlfn($phi, $s, cos($phi), $en) - $arg;
377
            //$phi -= $t * ($t * sqrt($t)) * $k;
378
            $t = ($this->pj_mlfn( $phi, $s, cos( $phi ), $en ) - $arg) * ($t * sqrt( $t )) * $k;
379
            $phi -= $t;
380
            if( abs( $t ) < Proj4php::$common->EPSLN )
381
                return $phi;
382
        }
383
 
384
        Proj4php::reportError( "cass:pj_inv_mlfn: Convergence error" );
385
 
386
        return $phi;
387
    }
388
 
389
}