SICOPOLIS V5-dev  Revision 1279
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-2018 Ralf Greve, Reinhard Calov, Alex Robinson, Fuyuki Saito
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  use sico_types_m
41  use error_m
42 
43  implicit none
44 
47 
48 contains
49 
50 !-------------------------------------------------------------------------------
51 !> Forward stereographic projection for an ellipsoidal planet.
52 !<------------------------------------------------------------------------------
53  subroutine stereo_forw_ellipsoid(lambda_val, phi_val, A, B, &
54  lambda0, phi0, x_val, y_val)
55 
56  implicit none
57 
58  real(dp), intent(in) :: lambda_val, phi_val, A, B, lambda0, phi0
59  real(dp), intent(out) :: x_val, y_val
60 
61  integer(i4b) :: l
62  integer(i4b) :: sign_phi0
63 
64  real(dp) :: phi_aux, phi0_aux
65  real(dp) :: e, mc, t, tc, kp, rho, phi_p
66  real(dp) :: sinphi0, sinlambda0, cosphi0, coslambda0
67 
68  if (phi0 > eps) then ! for northern hemisphere
69  sign_phi0 = 1
70  else if (phi0 < (-eps)) then ! for southern hemisphere
71  sign_phi0 = -1
72  else
73  errormsg = ' >>> stereo_forw_ellipsoid: phi0 must be different from zero!'
74  call error(errormsg)
75  end if
76 
77  phi_aux = phi_val * sign_phi0
78  phi0_aux = phi0 * sign_phi0
79 
80  e=sqrt((a**2-b**2)/(a**2))
81 
82  sinphi0 = sin(phi0_aux)
83  sinlambda0 = sin(lambda0)
84  cosphi0 = cos(phi0_aux)
85  coslambda0 = cos(lambda0)
86 
87  mc=cosphi0/sqrt(1.0_dp-e*e*sinphi0*sinphi0)
88  t=sqrt(((1.0_dp-sin(phi_aux))/(1.0_dp+sin(phi_aux)))* &
89  ((1.0_dp+e*sin(phi_aux))/(1.0_dp-e*sin(phi_aux)))**e)
90  tc=sqrt(((1.0_dp-sinphi0)/(1.0_dp+sinphi0))* &
91  ((1.0_dp+e*sinphi0)/(1.0_dp-e*sinphi0))**e)
92  rho=a*mc*t/tc
93 
94  x_val = rho*sin(lambda_val-lambda0)
95  y_val = -sign_phi0 * rho*cos(lambda_val-lambda0)
96 
97  end subroutine stereo_forw_ellipsoid
98 
99 !-------------------------------------------------------------------------------
100 !> Inverse stereographic projection for an ellipsoidal planet.
101 !<------------------------------------------------------------------------------
102  subroutine stereo_inv_ellipsoid(x_val, y_val, A, B, &
103  lambda0, phi0, lambda_val, phi_val)
105  implicit none
106 
107  real(dp), intent(in) :: x_val, y_val, A, B, lambda0, phi0
108  real(dp), intent(out) :: lambda_val, phi_val
109 
110  integer(i4b) :: l
111  integer(i4b) :: sign_phi0
112  real(dp) :: phi_aux, phi0_aux
113  real(dp) :: e, mc, t, tc, kp, rho, phi_p, residual
114  real(dp) :: sinphi0, sinlambda0, cosphi0, coslambda0
115 
116  real(dp), parameter :: eps_residual = 1.0e-09_dp
117 
118  if (phi0 > eps) then ! for northern hemisphere
119  sign_phi0 = 1
120  else if (phi0 < (-eps)) then ! for southern hemisphere
121  sign_phi0 = -1
122  else
123  errormsg = ' >>> stereo_inv_ellipsoid: phi0 must be different from zero!'
124  call error(errormsg)
125  end if
126 
127  phi0_aux = phi0 * sign_phi0
128 
129  e=sqrt((a**2-b**2)/(a**2))
130 
131  sinphi0 = sin(phi0_aux)
132  sinlambda0 = sin(lambda0)
133  cosphi0 = cos(phi0_aux)
134  coslambda0 = cos(lambda0)
135 
136  tc=sqrt(((1.0_dp-sinphi0)/(1.0_dp+sinphi0))* &
137  ((1.0_dp+e*sinphi0)/(1.0_dp-e*sinphi0))**e)
138  mc=cosphi0/sqrt(1.0_dp-e*e*sinphi0*sinphi0)
139  rho=sqrt(x_val*x_val+y_val*y_val)
140  t=rho*tc/(a*mc)
141 
142  if ((x_val /= 0.0_dp).or.(y_val /= 0.0_dp)) then
143  lambda_val = lambda0 + sign_phi0*atan2(y_val,x_val) + 0.5_dp*pi
144  else
145  lambda_val = lambda0 + 0.5_dp*pi
146  end if
147 
148  ! fix-point iteration
149 
150  phi_p=0.5_dp*pi-2.0_dp*atan(t)
151  l=0
152  residual=3600.0_dp
153  do while(residual >= eps_residual)
154  l=l+1
155  phi_aux=0.5_dp*pi-2.0_dp*atan(t*((1.0_dp-e*sin(phi_p))/ &
156  (1.0_dp+e*sin(phi_p)))**(0.5_dp*e))
157  residual=abs(phi_aux-phi_p)
158  phi_p=phi_aux
159  end do
160 
161  phi_val = phi_aux * sign_phi0
162 
163  if (lambda_val < 0.0_dp) then
164  lambda_val = lambda_val + 2.0_dp*pi
165  else if (lambda_val >= (2.0_dp*pi)) then
166  lambda_val = lambda_val - 2.0_dp*pi
167  end if
168 
169  end subroutine stereo_inv_ellipsoid
170 
171 !-------------------------------------------------------------------------------
172 !> Forward stereographic projection for a spherical planet.
173 !<------------------------------------------------------------------------------
174  subroutine stereo_forw_sphere(lambda_val, phi_val, R, lambda0, phi0, &
175  x_val, y_val)
177  implicit none
178 
179  real(dp), intent(in) :: lambda_val, phi_val, R, lambda0, phi0
180  real(dp), intent(out) :: x_val, y_val
181 
182  real(dp) :: K
183 
184  if (phi0 > eps) then ! for northern hemisphere
185 
186  k = (cos(0.25_dp*pi-0.5_dp*phi0))**2
187 
188  x_val = 2.0_dp*r*k*tan(0.25_dp*pi-0.5_dp*phi_val) &
189  *sin(lambda_val-lambda0)
190  y_val = -2.0_dp*r*k*tan(0.25_dp*pi-0.5_dp*phi_val) &
191  *cos(lambda_val-lambda0)
192 
193  else if (phi0 < (-eps)) then ! for southern hemisphere
194 
195  k = (cos(0.25_dp*pi+0.5_dp*phi0))**2
196 
197  x_val = 2.0_dp*r*k*tan(0.25_dp*pi+0.5_dp*phi_val) &
198  *sin(lambda_val-lambda0)
199  y_val = 2.0_dp*r*k*tan(0.25_dp*pi+0.5_dp*phi_val) &
200  *cos(lambda_val-lambda0)
201 
202  else
203  errormsg = ' >>> stereo_forw_sphere: phi0 must be different from zero!'
204  call error(errormsg)
205  end if
206 
207  end subroutine stereo_forw_sphere
208 
209 !-------------------------------------------------------------------------------
210 !> Inverse stereographic projection for a spherical planet.
211 !<------------------------------------------------------------------------------
212  subroutine stereo_inv_sphere(x_val, y_val, R, lambda0, phi0, &
213  lambda_val, phi_val)
215  implicit none
216 
217  real(dp), intent(in) :: x_val, y_val, R, lambda0, phi0
218  real(dp), intent(out) :: lambda_val, phi_val
219 
220  real(dp) :: K
221 
222  if (phi0 > eps) then ! for northern hemisphere
223  k = (cos(0.25_dp*pi-0.5_dp*phi0))**2
224 
225  phi_val = 0.5_dp*pi &
226  -2.0_dp*atan(sqrt(x_val**2+y_val**2)/(2.0_dp*r*k))
227 
228  if ((x_val /= 0.0_dp).or.(y_val /= 0.0_dp)) then
229  lambda_val = lambda0 + atan2(y_val,x_val) + 0.5_dp*pi
230  else
231  lambda_val = lambda0 + 0.5_dp*pi
232  end if
233 
234  if (lambda_val < 0.0_dp) then
235  lambda_val = lambda_val + 2.0_dp*pi
236  else if (lambda_val >= (2.0_dp*pi)) then
237  lambda_val = lambda_val - 2.0_dp*pi
238  end if
239 
240  else if (phi0 < (-eps)) then ! for southern hemisphere
241 
242  k = (cos(0.25_dp*pi+0.5_dp*phi0))**2
243 
244  phi_val = -0.5_dp*pi &
245  +2.0_dp*atan(sqrt(x_val**2+y_val**2)/(2.0_dp*r*k))
246 
247  if ((x_val /= 0.0_dp).or.(y_val /= 0.0_dp)) then
248  lambda_val = lambda0 - atan2(y_val,x_val) + 0.5_dp*pi
249  else
250  lambda_val = lambda0 + 0.5_dp*pi
251  end if
252 
253  if (lambda_val < 0.0_dp) then
254  lambda_val = lambda_val + 2.0_dp*pi
255  else if (lambda_val >= (2.0_dp*pi)) then
256  lambda_val = lambda_val - 2.0_dp*pi
257  end if
258  else
259  errormsg = ' >>> stereo_inv_sphere: phi0 must be different from zero!'
260  call error(errormsg)
261  end if
262 
263  end subroutine stereo_inv_sphere
264 
265 !-------------------------------------------------------------------------------
266 
267 end module stereo_proj_m
268 !
subroutine, public stereo_inv_sphere(x_val, y_val, R, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for a spherical planet.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
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.
subroutine, public stereo_forw_ellipsoid(lambda_val, phi_val, A, B, lambda0, phi0, x_val, y_val)
Forward stereographic projection for an ellipsoidal planet.
Writing of error messages and stopping execution.
Definition: error_m.F90:35
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
real(dp), parameter pi
pi: Constant pi
character(len=256) errormsg
errormsg: Error-message string
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
Declarations of global variables for SICOPOLIS.