SICOPOLIS V5-dev  Revision 1173
stereo_proj_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s t e r e o _ p r o j _ m
4 !
5 !> @file
6 !!
7 !! Computation of the forward or inverse stereographic projection,
8 !! alternatively for a spherical or an ellipsoidal planet.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2017 Ralf Greve, Reinhard Calov, Alex Robinson
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Computation of the forward or inverse stereographic projection,
35 !! alternatively for a spherical or an ellipsoidal planet.
36 !<------------------------------------------------------------------------------
38 
39 #if (defined(MODEL_SICOPOLIS))
40  use sico_types_m, only : i4b, dp
41  use sico_variables_m, only : pi, eps
42 #endif
43 
44  implicit none
45 
46  private
49 
50 #if (!defined(MODEL_SICOPOLIS))
51  integer, private, parameter :: i4b = selected_int_kind(9) ! 4-byte integers
52  integer, private, parameter :: dp = kind(1.0d0) ! double-precision reals
53  real(dp), private, parameter :: pi = 3.141592653589793_dp
54  real(dp), private, parameter :: eps = 1.0e-05_dp
55 #endif
56 
57 contains
58 
59 !-------------------------------------------------------------------------------
60 !> Forward stereographic projection for an ellipsoidal planet.
61 !<------------------------------------------------------------------------------
62  subroutine stereo_forw_ellipsoid(lambda_val, phi_val, A, B, &
63  lambda0, phi0, x_val, y_val)
64 
65  implicit none
66 
67  real(dp), intent(in) :: lambda_val, phi_val, A, B, lambda0, phi0
68  real(dp), intent(out) :: x_val, y_val
69 
70  integer(i4b) :: l
71  integer(i4b) :: sign_phi0
72  real(dp) :: phi_aux, phi0_aux
73  real(dp) :: e, mc, t, tc, kp, rho, phi_p
74  real(dp) :: sinphi0, sinlambda0, cosphi0, coslambda0
75 
76  if (phi0 > eps) then ! for northern hemisphere
77  sign_phi0 = 1
78  else if (phi0 < (-eps)) then ! for southern hemisphere
79  sign_phi0 = -1
80  else
81  stop ' >>> stereo_forw_ellipsoid: phi0 must be different from zero!'
82  end if
83 
84  phi_aux = phi_val * sign_phi0
85  phi0_aux = phi0 * sign_phi0
86 
87  e=sqrt((a**2-b**2)/(a**2))
88 
89  sinphi0 = sin(phi0_aux)
90  sinlambda0 = sin(lambda0)
91  cosphi0 = cos(phi0_aux)
92  coslambda0 = cos(lambda0)
93 
94  mc=cosphi0/sqrt(1.0_dp-e*e*sinphi0*sinphi0)
95  t=sqrt(((1.0_dp-sin(phi_aux))/(1.0_dp+sin(phi_aux)))* &
96  ((1.0_dp+e*sin(phi_aux))/(1.0_dp-e*sin(phi_aux)))**e)
97  tc=sqrt(((1.0_dp-sinphi0)/(1.0_dp+sinphi0))* &
98  ((1.0_dp+e*sinphi0)/(1.0_dp-e*sinphi0))**e)
99  rho=a*mc*t/tc
100 
101  x_val = rho*sin(lambda_val-lambda0)
102  y_val = -sign_phi0 * rho*cos(lambda_val-lambda0)
103 
104  end subroutine stereo_forw_ellipsoid
105 
106 !-------------------------------------------------------------------------------
107 !> Inverse stereographic projection for an ellipsoidal planet.
108 !<------------------------------------------------------------------------------
109  subroutine stereo_inv_ellipsoid(x_val, y_val, A, B, &
110  lambda0, phi0, lambda_val, phi_val)
112  implicit none
113 
114  real(dp), intent(in) :: x_val, y_val, A, B, lambda0, phi0
115  real(dp), intent(out) :: lambda_val, phi_val
116 
117  integer(i4b) :: l
118  integer(i4b) :: sign_phi0
119  real(dp) :: phi_aux, phi0_aux
120  real(dp) :: e, mc, t, tc, kp, rho, phi_p, residual
121  real(dp) :: sinphi0, sinlambda0, cosphi0, coslambda0
122 
123  real(dp), parameter :: eps_residual = 1.0e-09_dp
124 
125  if (phi0 > eps) then ! for northern hemisphere
126  sign_phi0 = 1
127  else if (phi0 < (-eps)) then ! for southern hemisphere
128  sign_phi0 = -1
129  else
130  stop ' >>> stereo_inv_ellipsoid: phi0 must be different from zero!'
131  end if
132 
133  phi0_aux = phi0 * sign_phi0
134 
135  e=sqrt((a**2-b**2)/(a**2))
136 
137  sinphi0 = sin(phi0_aux)
138  sinlambda0 = sin(lambda0)
139  cosphi0 = cos(phi0_aux)
140  coslambda0 = cos(lambda0)
141 
142  tc=sqrt(((1.0_dp-sinphi0)/(1.0_dp+sinphi0))* &
143  ((1.0_dp+e*sinphi0)/(1.0_dp-e*sinphi0))**e)
144  mc=cosphi0/sqrt(1.0_dp-e*e*sinphi0*sinphi0)
145  rho=sqrt(x_val*x_val+y_val*y_val)
146  t=rho*tc/(a*mc)
147 
148  lambda_val = lambda0 + sign_phi0*atan2(y_val,x_val) + 0.5_dp*pi
149 
150  ! fix-point iteration
151 
152  phi_p=0.5_dp*pi-2.0_dp*atan(t)
153  l=0
154  residual=3600.0_dp
155  do while(residual >= eps_residual)
156  l=l+1
157  phi_aux=0.5_dp*pi-2.0_dp*atan(t*((1.0_dp-e*sin(phi_p))/ &
158  (1.0_dp+e*sin(phi_p)))**(0.5_dp*e))
159  residual=abs(phi_aux-phi_p)
160  phi_p=phi_aux
161  end do
162 
163  phi_val = phi_aux * sign_phi0
164 
165  if (lambda_val < 0.0_dp) then
166  lambda_val = lambda_val + 2.0_dp*pi
167  else if (lambda_val >= (2.0_dp*pi)) then
168  lambda_val = lambda_val - 2.0_dp*pi
169  end if
170 
171  end subroutine stereo_inv_ellipsoid
172 
173 !-------------------------------------------------------------------------------
174 !> Forward stereographic projection for a spherical planet.
175 !<------------------------------------------------------------------------------
176  subroutine stereo_forw_sphere(lambda_val, phi_val, R, lambda0, phi0, &
177  x_val, y_val)
179  implicit none
180 
181  real(dp), intent(in) :: lambda_val, phi_val, R, lambda0, phi0
182  real(dp), intent(out) :: x_val, y_val
183 
184  real(dp) :: K
185 
186  if (phi0 > eps) then ! for northern hemisphere
187 
188  k = (cos(0.25_dp*pi-0.5_dp*phi0))**2
189 
190  x_val = 2.0_dp*r*k*tan(0.25_dp*pi-0.5_dp*phi_val) &
191  *sin(lambda_val-lambda0)
192  y_val = -2.0_dp*r*k*tan(0.25_dp*pi-0.5_dp*phi_val) &
193  *cos(lambda_val-lambda0)
194 
195  else if (phi0 < (-eps)) then ! for southern hemisphere
196 
197  k = (cos(0.25_dp*pi+0.5_dp*phi0))**2
198 
199  x_val = 2.0_dp*r*k*tan(0.25_dp*pi+0.5_dp*phi_val) &
200  *sin(lambda_val-lambda0)
201  y_val = 2.0_dp*r*k*tan(0.25_dp*pi+0.5_dp*phi_val) &
202  *cos(lambda_val-lambda0)
203 
204  else
205 
206  stop ' >>> stereo_forw_sphere: phi0 must be different from zero!'
207 
208  end if
209 
210  end subroutine stereo_forw_sphere
211 
212 !-------------------------------------------------------------------------------
213 !> Inverse stereographic projection for a spherical planet.
214 !<------------------------------------------------------------------------------
215  subroutine stereo_inv_sphere(x_val, y_val, R, lambda0, phi0, &
216  lambda_val, phi_val)
218  implicit none
219 
220  real(dp), intent(in) :: x_val, y_val, R, lambda0, phi0
221  real(dp), intent(out) :: lambda_val, phi_val
222 
223  real(dp) :: K
224 
225  if (phi0 > eps) then ! for northern hemisphere
226 
227  k = (cos(0.25_dp*pi-0.5_dp*phi0))**2
228 
229  phi_val = 0.5_dp*pi &
230  -2.0_dp*atan(sqrt(x_val**2+y_val**2)/(2.0_dp*r*k))
231  lambda_val = lambda0 + atan2(y_val,x_val)+0.5_dp*pi
232 
233  if (lambda_val < 0.0_dp) then
234  lambda_val = lambda_val + 2.0_dp*pi
235  else if (lambda_val >= (2.0_dp*pi)) then
236  lambda_val = lambda_val - 2.0_dp*pi
237  end if
238 
239  else if (phi0 < (-eps)) then ! for southern hemisphere
240 
241  k = (cos(0.25_dp*pi+0.5_dp*phi0))**2
242 
243  phi_val = -0.5_dp*pi &
244  +2.0_dp*atan(sqrt(x_val**2+y_val**2)/(2.0_dp*r*k))
245  lambda_val = lambda0 - atan2(y_val,x_val)+0.5_dp*pi
246 
247  if (lambda_val < 0.0_dp) then
248  lambda_val = lambda_val + 2.0_dp*pi
249  else if (lambda_val >= (2.0_dp*pi)) then
250  lambda_val = lambda_val - 2.0_dp*pi
251  end if
252 
253  else
254 
255  stop ' >>> stereo_inv_sphere: phi0 must be different from zero!'
256 
257  end if
258 
259  end subroutine stereo_inv_sphere
260 
261 !-------------------------------------------------------------------------------
262 
263 end module stereo_proj_m
264 !
subroutine, public stereo_inv_sphere(x_val, y_val, R, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for a spherical planet.
subroutine, public stereo_forw_sphere(lambda_val, phi_val, R, lambda0, phi0, x_val, y_val)
Forward stereographic projection for a spherical planet.
real(dp), parameter eps
eps: Small number
subroutine, public stereo_inv_ellipsoid(x_val, y_val, A, B, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for an ellipsoidal planet.
Declarations of kind types for SICOPOLIS.
integer, parameter i4b
4-byte integers
subroutine, public stereo_forw_ellipsoid(lambda_val, phi_val, A, B, lambda0, phi0, x_val, y_val)
Forward stereographic projection for an ellipsoidal planet.
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
real(dp), parameter pi
pi: Constant pi
integer, parameter dp
Double-precision reals.
Declarations of global variables for SICOPOLIS.