Subversion Repositories eFlore/Projets.eflore-projets

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
1001 delphine 1
<?php
2
/**
3
 * Author : Julien Moquet
4
 *
5
 * Inspired by Proj4php from Mike Adair madairATdmsolutions.ca
6
 *                      and Richard Greenwood rich@greenwoodma$p->com
7
 * License: LGPL as per: http://www.gnu.org/copyleft/lesser.html
8
 */
9
 
10
// Initialize the Stereographic projection
11
class Proj4phpProjStere {
12
 
13
    protected $TOL = 1.e-8;
14
    protected $NITER = 8;
15
    protected $CONV = 1.e-10;
16
    protected $S_POLE = 0;
17
    protected $N_POLE = 1;
18
    protected $OBLIQ = 2;
19
    protected $EQUIT = 3;
20
 
21
    /**
22
     *
23
     * @param type $phit
24
     * @param type $sinphi
25
     * @param type $eccen
26
     * @return type
27
     */
28
    public function ssfn_( $phit, $sinphi, $eccen ) {
29
        $sinphi *= $eccen;
30
        return (tan( .5 * (Proj4php::$common->HALF_PI + $phit) ) * pow( (1. - $sinphi) / (1. + $sinphi), .5 * $eccen ));
31
    }
32
 
33
    /**
34
     *
35
     */
36
    public function init() {
37
        $this->phits = $this->lat_ts ? $this->lat_ts : Proj4php::$common->HALF_PI;
38
        $t = abs( $this->lat0 );
39
        if( (abs( $t ) - Proj4php::$common->HALF_PI) < Proj4php::$common->EPSLN ) {
40
            $this->mode = $this->lat0 < 0. ? $this->S_POLE : $this->N_POLE;
41
        } else {
42
            $this->mode = $t > Proj4php::$common->EPSLN ? $this->OBLIQ : $this->EQUIT;
43
        }
44
        $this->phits = abs( $this->phits );
45
        if( $this->es ) {
46
            #$X;
47
 
48
            switch( $this->mode ) {
49
                case $this->N_POLE:
50
                case $this->S_POLE:
51
                    if( abs( $this->phits - Proj4php::$common->HALF_PI ) < Proj4php::$common->EPSLN ) {
52
                        $this->akm1 = 2. * $this->k0 / sqrt( pow( 1 + $this->e, 1 + $this->e ) * pow( 1 - $this->e, 1 - $this->e ) );
53
                    } else {
54
                        $t = sin( $this->phits );
55
                        $this->akm1 = cos( $this->phits ) / Proj4php::$common->tsfnz( $this->e, $this->phits, $t );
56
                        $t *= $this->e;
57
                        $this->akm1 /= sqrt( 1. - $t * $t );
58
                    }
59
                    break;
60
                case $this->EQUIT:
61
                    $this->akm1 = 2. * $this->k0;
62
                    break;
63
                case $this->OBLIQ:
64
                    $t = sin( $this->lat0 );
65
                    $X = 2. * atan( $this->ssfn_( $this->lat0, $t, $this->e ) ) - Proj4php::$common->HALF_PI;
66
                    $t *= $this->e;
67
                    $this->akm1 = 2. * $this->k0 * cos( $this->lat0 ) / sqrt( 1. - $t * $t );
68
                    $this->sinX1 = sin( $X );
69
                    $this->cosX1 = cos( $X );
70
                    break;
71
            }
72
        } else {
73
            switch( $this->mode ) {
74
                case $this->OBLIQ:
75
                    $this->sinph0 = sin( $this->lat0 );
76
                    $this->cosph0 = cos( $this->lat0 );
77
                case $this->EQUIT:
78
                    $this->akm1 = 2. * $this->k0;
79
                    break;
80
                case $this->S_POLE:
81
                case $this->N_POLE:
82
                    $this->akm1 = abs( $this->phits - Proj4php::$common->HALF_PI ) >= Proj4php::$common->EPSLN ?
83
                              cos( $this->phits ) / tan( Proj4php::$common->FORTPI - .5 * $this->phits ) :
84
                              2. * $this->k0;
85
                    break;
86
            }
87
        }
88
    }
89
 
90
    /**
91
     * Stereographic forward equations--mapping lat,long to x,y
92
     *
93
     * @param type $p
94
     * @return type
95
     */
96
    public function forward( $p ) {
97
 
98
        $lon = $p->x;
99
        $lon = Proj4php::$common->adjust_lon( $lon - $this->long0 );
100
        $lat = $p->y;
101
        #$x;
102
        #$y;
103
 
104
        if( $this->sphere ) {
105
            /*
106
            $sinphi;
107
            $cosphi;
108
            $coslam;
109
            $sinlam;
110
            */
111
 
112
            $sinphi = sin( $lat );
113
            $cosphi = cos( $lat );
114
            $coslam = cos( $lon );
115
            $sinlam = sin( $lon );
116
            switch( $this->mode ) {
117
                case $this->EQUIT:
118
                    $y = 1. + $cosphi * $coslam;
119
                    if( y <= Proj4php::$common->EPSLN ) {
120
                        Proj4php::reportError("stere:forward:Equit");
121
                    }
122
                    $y = $this->akm1 / $y;
123
                    $x = $y * $cosphi * $sinlam;
124
                    $y *= $sinphi;
125
                    break;
126
                case $this->OBLIQ:
127
                    $y = 1. + $this->sinph0 * $sinphi + $this->cosph0 * $cosphi * $coslam;
128
                    if( $y <= Proj4php::$common->EPSLN ) {
129
                        Proj4php::reportError("stere:forward:Obliq");
130
                    }
131
                    $y = $this->akm1 / $y;
132
                    $x = $y * $cosphi * $sinlam;
133
                    $y *= $this->cosph0 * $sinphi - $this->sinph0 * $cosphi * $coslam;
134
                    break;
135
                case $this->N_POLE:
136
                    $coslam = -$coslam;
137
                    $lat = -$lat;
138
                //Note  no break here so it conitnues through S_POLE
139
                case $this->S_POLE:
140
                    if( abs( $lat - Proj4php::$common->HALF_PI ) < $this->TOL ) {
141
                        Proj4php::reportError("stere:forward:S_POLE");
142
                    }
143
                    $y = $this->akm1 * tan( Proj4php::$common->FORTPI + .5 * $lat );
144
                    $x = $sinlam * $y;
145
                    $y *= $coslam;
146
                    break;
147
            }
148
        } else {
149
            $coslam = cos( $lon );
150
            $sinlam = sin( $lon );
151
            $sinphi = sin( $lat );
152
            if( $this->mode == $this->OBLIQ || $this->mode == $this->EQUIT ) {
153
                $Xt = 2. * atan( $this->ssfn_( $lat, $sinphi, $this->e ) );
154
                $sinX = sin( $Xt - Proj4php::$common->HALF_PI );
155
                $cosX = cos( $Xt );
156
            }
157
            switch( $this->mode ) {
158
                case $this->OBLIQ:
159
                    $A = $this->akm1 / ($this->cosX1 * (1. + $this->sinX1 * $sinX + $this->cosX1 * $cosX * $coslam));
160
                    $y = $A * ($this->cosX1 * $sinX - $this->sinX1 * $cosX * $coslam);
161
                    $x = $A * $cosX;
162
                    break;
163
                case $this->EQUIT:
164
                    $A = 2. * $this->akm1 / (1. + $cosX * $coslam);
165
                    $y = $A * $sinX;
166
                    $x = $A * $cosX;
167
                    break;
168
                case $this->S_POLE:
169
                    $lat = -$lat;
170
                    $coslam = - $coslam;
171
                    $sinphi = -$sinphi;
172
                case $this->N_POLE:
173
                    $x = $this->akm1 * Proj4php::$common->tsfnz( $this->e, $lat, $sinphi );
174
                    $y = - $x * $coslam;
175
                    break;
176
            }
177
            $x = $x * $sinlam;
178
        }
179
 
180
        $p->x = $x * $this->a + $this->x0;
181
        $p->y = $y * $this->a + $this->y0;
182
 
183
        return $p;
184
    }
185
 
186
 
187
    /**
188
     * Stereographic inverse equations--mapping x,y to lat/long
189
     *
190
     * @param type $p
191
     * @return type
192
     */
193
    public function inverse( $p ) {
194
        $x = ($p->x - $this->x0) / $this->a;   /* descale and de-offset */
195
        $y = ($p->y - $this->y0) / $this->a;
196
        /*
197
        $lon;
198
        $lat;
199
        $cosphi;
200
        $sinphi;
201
        $rho;
202
        $tp = 0.0;
203
        $phi_l = 0.0;
204
        $i;
205
        */
206
        $halfe = 0.0;
207
        $pi2 = 0.0;
208
 
209
        if( $this->sphere ) {
210
            /*
211
            $c;
212
            $rh;
213
            $sinc;
214
            $cosc;
215
            */
216
 
217
            $rh = sqrt( $x * $x + $y * $y );
218
            $c = 2. * atan( $rh / $this->akm1 );
219
            $sinc = sin( $c );
220
            $cosc = cos( $c );
221
            $lon = 0.;
222
 
223
            switch( $this->mode ) {
224
                case $this->EQUIT:
225
                    if( abs( $rh ) <= Proj4php::$common->EPSLN ) {
226
                        $lat = 0.;
227
                    } else {
228
                        $lat = asin( $y * $sinc / $rh );
229
                    }
230
                    if( $cosc != 0. || $x != 0. )
231
                        $lon = atan2( $x * $sinc, $cosc * $rh );
232
                    break;
233
                case $this->OBLIQ:
234
                    if( abs( $rh ) <= Proj4php::$common->EPSLN ) {
235
                        $lat = $this->phi0;
236
                    } else {
237
                        $lat = asin( $cosc * $this->sinph0 + $y * $sinc * $this->cosph0 / $rh );
238
                    }
239
                    $c = $cosc - $this->sinph0 * sin( $lat );
240
                    if( $c != 0. || $x != 0. ) {
241
                        $lon = atan2( $x * $sinc * $this->cosph0, $c * $rh );
242
                    }
243
                    break;
244
                case $this->N_POLE:
245
                    $y = -$y;
246
                case $this->S_POLE:
247
                    if( abs( $rh ) <= Proj4php::$common->EPSLN ) {
248
                        $lat = $this->phi0;
249
                    } else {
250
                        $lat = asin( $this->mode == $this->S_POLE ? -$cosc : $cosc  );
251
                    }
252
                    $lon = ($x == 0. && $y == 0.) ? 0. : atan2( $x, $y );
253
                    break;
254
            }
255
            $p->x = Proj4php::$common->adjust_lon( $lon + $this->long0 );
256
            $p->y = $lat;
257
        } else {
258
            $rho = sqrt( $x * $x + $y * $y );
259
            switch( $this->mode ) {
260
                case $this->OBLIQ:
261
                case $this->EQUIT:
262
                    $tp = 2. * atan2( $rho * $this->cosX1, $this->akm1 );
263
                    $cosphi = cos( $tp );
264
                    $sinphi = sin( $tp );
265
                    if( $rho == 0.0 ) {
266
                        $phi_l = asin( $cosphi * $this->sinX1 );
267
                    } else {
268
                        $phi_l = asin( $cosphi * $this->sinX1 + ($y * $sinphi * $this->cosX1 / $rho) );
269
                    }
270
 
271
                    $tp = tan( .5 * (Proj4php::$common->HALF_PI + $phi_l) );
272
                    $x *= $sinphi;
273
                    $y = $rho * $this->cosX1 * $cosphi - $y * $this->sinX1 * $sinphi;
274
                    $pi2 = Proj4php::$common->HALF_PI;
275
                    $halfe = .5 * $this->e;
276
                    break;
277
                case $this->N_POLE:
278
                    $y = -$y;
279
                case $this->S_POLE:
280
                    $tp = - $rho / $this->akm1;
281
                    $phi_l = Proj4php::$common->HALF_PI - 2. * atan( $tp );
282
                    $pi2 = -Proj4php::$common->HALF_PI;
283
                    $halfe = -.5 * $this->e;
284
                    break;
285
            }
286
            for( $i = $this->NITER; $i--; $phi_l = $lat ) { //check this
287
                $sinphi = $this->e * sin( $phi_l );
288
                $lat = 2. * atan( $tp * pow( (1. + $sinphi) / (1. - $sinphi), $halfe ) ) - $pi2;
289
                if( abs( phi_l - lat ) < $this->CONV ) {
290
                    if( $this->mode == $this->S_POLE )
291
                        $lat = -$lat;
292
                    $lon = ($x == 0. && $y == 0.) ? 0. : atan2( $x, $y );
293
                    $p->x = Proj4php::$common->adjust_lon( $lon + $this->long0 );
294
                    $p->y = $lat;
295
                    return $p;
296
                }
297
            }
298
        }
299
    }
300
 
301
}
302
 
303
Proj4php::$proj['stere'] = new Proj4phpProjStere();