Ticket #325: azi.f90

File azi.f90, 9.6 KB (added by Ian Culverwell, 11 years ago)
Line 
1!ifort -o azi.exe azi.f90
2!
3PROGRAM azi
4
5 IMPLICIT none
6
7 INTEGER, PARAMETER :: wp=KIND(1.0D0)
8 REAL(KIND=wp), PARAMETER :: pi=4.d0*ATAN(1.d0), dtor=pi/180.d0, rad2deg=1.0d0/dtor
9 REAL(KIND=wp), PARAMETER :: r_earth=6.4d6
10
11 INTEGER, PARAMETER :: npts=1
12 INTEGER :: i
13
14 REAL(KIND=wp), DIMENSION(npts,3) :: r_leo ! in ECF coords
15 REAL(KIND=wp), DIMENSION(npts) :: rad_leo, lon_leo, lat_leo ! in ECF coords
16
17 REAL(KIND=wp), DIMENSION(npts,3) :: r_gns ! in ECF coords
18 REAL(KIND=wp), DIMENSION(npts) :: rad_gns, lon_gns, lat_gns ! in ECF coords
19
20 REAL(KIND=wp), DIMENSION(3), PARAMETER :: pa=(/ 0.0d0, 0.0d0, 1.0d0 /)
21 REAL(KIND=wp), DIMENSION(npts,3) :: perigee
22 REAL(KIND=wp) :: slta, ro, alpha, theta
23 REAL(KIND=wp), DIMENSION(npts) :: azimuth_tp, azimuth_ropp, azimuth_gpac
24
25 CHARACTER(LEN=256) :: title
26
27! 1.1 Define satellite positions (ECF coordinates)
28! ------------------------------------------------
29
30 DO i=1,npts
31
32! title = 'LEO due west of GNS'
33! rad_leo(i) = r_earth + 800.0d3 ; lat_leo(i) = 0.0d0 * dtor ; lon_leo(i) = -10.d0 * dtor
34! rad_gns(i) = r_earth + 20000.0d3 ; lat_gns(i) = 0.0d0 * dtor ; lon_gns(i) = 90.d0 * dtor
35
36! title = 'LEO due east of GNS'
37! rad_leo(i) = r_earth + 800.0d3 ; lat_leo(i) = 0.0d0 * dtor ; lon_leo(i) = 10.d0 * dtor
38! rad_gns(i) = r_earth + 20000.0d3 ; lat_gns(i) = 0.0d0 * dtor ; lon_gns(i) = -90.d0 * dtor
39
40 title = 'LEO due north of GNS'
41 rad_leo(i) = r_earth + 800.0d3 ; lat_leo(i) = 10.0d0 * dtor ; lon_leo(i) = 0.d0 * dtor
42 rad_gns(i) = r_earth + 20000.0d3 ; lat_gns(i) = -10.0d0 * dtor ; lon_gns(i) = 0.d0 * dtor
43
44! title = 'LEO due south of GNS'
45! rad_leo(i) = r_earth + 800.0d3 ; lat_leo(i) = -10.0d0 * dtor ; lon_leo(i) = 0.d0 * dtor
46! rad_gns(i) = r_earth + 20000.0d3 ; lat_gns(i) = 10.0d0 * dtor ; lon_gns(i) = 0.d0 * dtor
47
48! title = 'LEO and GEO as in 1st picture of #233'
49! rad_leo(i) = r_earth + 800.0d3 ; lat_leo(i) = 55.0d0 * dtor ; lon_leo(i) = -25.d0 * dtor
50! rad_gns(i) = r_earth + 20000.0d3 ; lat_gns(i) = -50.0d0 * dtor ; lon_gns(i) = 0.d0 * dtor
51
52! title = 'LEO and GEO as in 2nd picture of #233'
53! rad_leo(i) = r_earth + 800.0d3 ; lat_leo(i) = 60.0d0 * dtor ; lon_leo(i) = -40.d0 * dtor
54! rad_gns(i) = r_earth + 20000.0d3 ; lat_gns(i) = 5.0d0 * dtor ; lon_gns(i) = -170.d0 * dtor
55
56
57 r_leo(i,:) = rad_leo(i) * (/ COS(lon_leo(i))*COS(lat_leo(i)), SIN(lon_leo(i))*COS(lat_leo(i)), SIN(lat_leo(i)) /)
58
59 r_gns(i,:) = rad_gns(i) * (/ COS(lon_gns(i))*COS(lat_gns(i)), SIN(lon_gns(i))*COS(lat_gns(i)), SIN(lat_gns(i)) /)
60
61
62! 1.2 Determine ray tangent points
63! --------------------------------
64
65 slta = impact_parameter(r_leo(i,:), r_gns(i,:))
66 ro = SQRT(DOT_PRODUCT(r_leo(i,:), r_leo(i,:)))
67 alpha = ACOS(slta/ro)
68
69 perigee(i,:) = rotate(r_leo(i,:), vector_product(r_leo(i,:), r_gns(i,:)), alpha) * (slta/ro)
70
71 PRINT*,'slta,ro,alpha,perigee = ', slta,ro,alpha,perigee(i,:)
72
73 ENDDO
74
75
76! 1.4 Cross-section azimuth at tangent point as in ROPP
77! -----------------------------------------------------
78
79 DO i=1,size(r_leo,1)
80
81 theta = vector_angle(vector_product(r_gns(i,:),r_leo(i,:)), &
82 vector_product(pa, perigee(i,:)))
83
84 PRINT*,'theta_ropp = ', theta
85
86 PRINT*,'DOT_PRODUCT(r_gns(i,:)-r_leo(i,:),vector_product(pa, perigee(i,:))) = ', &
87 DOT_PRODUCT(r_gns(i,:)-r_leo(i,:),vector_product(pa, perigee(i,:)))
88
89 if (DOT_PRODUCT(r_gns(i,:)-r_leo(i,:),vector_product(pa, perigee(i,:))) < 0) THEN
90 theta = 2.0d0*Pi - theta
91 endif
92
93 PRINT*,'theta_ropp = ', theta
94
95 azimuth_tp(i) = theta * rad2deg
96
97 azimuth_ropp(i) = azimuth_tp(i)
98
99 ENDDO
100
101! 1.4 Cross-section azimuth at tangent point as in GPAC
102! -----------------------------------------------------
103
104 DO i=1,size(r_leo,1)
105
106 theta = vector_angle(vector_product(perigee(i,:), pa), &
107 vector_product(r_gns(i,:),r_leo(i,:)), -perigee(i,:))
108
109 PRINT*,'theta_gpac = ', theta
110
111 azimuth_tp(i) = theta * rad2deg
112
113 PRINT*,'azimuth_gpac = ', azimuth_tp(i)
114
115 if (azimuth_tp(i) < 0.0_wp ) azimuth_tp(i) = azimuth_tp(i) + 360.0_wp
116
117 PRINT*,'azimuth_gpac = ', azimuth_tp(i)
118
119 azimuth_gpac(i) = azimuth_tp(i)
120
121 ENDDO
122
123
124PRINT*, '*** ' // TRIM(title) // ' ***'
125PRINT*, 'azimuth_ropp = ', azimuth_ropp
126PRINT*, 'azimuth_gpac = ', azimuth_gpac
127
128
129
130CONTAINS
131
132!****f* Coordinates/rotate
133!
134! NAME
135! rotate - Rotate a vector in cartesian coordinates around
136! a given axis by a given angle
137!
138! SYNOPSIS
139! Rotate = rotate(X, A, phi)
140!
141! DESCRIPTION
142! This function rotates a vector X around a given axis A by angle phi.
143! N*(N,X) + [N,X]*Sin(Phi) + (X-N*(N,X))*Cos(Phi), where N=A/|A|
144!
145! INPUTS
146! X Vector to rotate
147! A Rotation axis
148! Phi Rotation angle (rad)
149!
150! OUTPUT
151! Rotate Rotated vector
152!
153! NOTES
154!
155! SEE ALSO
156!
157! REFERENCES
158!
159! AUTHOR
160! Met Office, Exeter, UK.
161! Any comments on this software should be given via the ROM SAF
162! Helpdesk at http://www.romsaf.org
163!
164! COPYRIGHT
165! (c) EUMETSAT. All rights reserved.
166! For further details please refer to the file COPYRIGHT
167! which you should have received as part of this distribution.
168!
169!****
170
171function rotate(X, A, Phi) result(R)
172
173! 1.1 Declarations
174! ----------------
175
176 implicit none
177
178 real(wp), dimension(3), intent(in) :: X ! input vector
179 real(wp), dimension(3), intent(in) :: A ! rotation axis
180 real(wp), intent(in) :: phi ! rotation angle
181 real(wp), dimension(3) :: R ! rotated vector
182
183 real(wp), dimension(3) :: norm ! normed rotation axis
184
185! 1.2 Frame rotation
186! ------------------
187
188! N*(N,X) + [N,X]*Sin(Phi) + (X-N*(N,X))*Cos(Phi), where N=A/|A|
189
190 norm = A(:)/Sqrt(Dot_Product(A(:), A(:)))
191
192 R = norm*(Dot_Product(norm, X)) + vector_product(norm, X)*sin(phi) &
193 + (X - norm*Dot_Product(norm,X))*cos(phi)
194
195end function rotate
196
197!****f* Coordinates/vector_product
198!
199! NAME
200! vector_product - Compute a vector product of two cartesian vectors
201!
202! SYNOPSIS
203! product = vector_product(X, Y)
204!
205! INPUTS
206! X Vector 1
207! Y Vector 2
208!
209! OUTPUT
210! Product Vector product
211!
212! AUTHOR
213! Met Office, Exeter, UK.
214! Any comments on this software should be given via the ROM SAF
215! Helpdesk at http://www.romsaf.org
216!
217! COPYRIGHT
218! (c) EUMETSAT. All rights reserved.
219! For further details please refer to the file COPYRIGHT
220! which you should have received as part of this distribution.
221!
222!****
223
224 function vector_product(X, Y) result(product)
225
226 real(wp), dimension(3), intent(in) :: X
227 real(wp), dimension(3), intent(in) :: Y
228 real(wp), dimension(3) :: product
229
230 product = (/ X(2)*Y(3) - X(3)*Y(2), &
231 X(3)*Y(1) - X(1)*Y(3), &
232 X(1)*Y(2) - X(2)*Y(1) /)
233
234 end function vector_product
235
236!****f* Coordinates/vector_angle
237!
238! NAME
239! vector_angle - Find the angle between two cartesian vectors
240!
241! SYNOPSIS
242! angle = vector_angle(X, Y, A)
243!
244! INPUTS
245! X Vector 1
246! Y Vector 2
247! A Orientation axis (optional)
248!
249! OUTPUT
250! Angle Angle between vectors
251!
252! AUTHOR
253! Met Office, Exeter, UK.
254! Any comments on this software should be given via the ROM SAF
255! Helpdesk at http://www.romsaf.org
256!
257! COPYRIGHT
258! (c) EUMETSAT. All rights reserved.
259! For further details please refer to the file COPYRIGHT
260! which you should have received as part of this distribution.
261!
262!****
263
264function vector_angle(X, Y, A) result(angle)
265
266 real(wp), dimension(3), intent(in) :: X
267 real(wp), dimension(3), intent(in) :: Y
268 real(wp), dimension(3), optional, intent(in) :: A
269 real(wp) :: angle
270
271 real(wp), dimension(3) :: n, alpha, beta, gamma
272 real(wp) :: nn
273
274 if (present(A)) then
275 n = A
276 else
277 n = vector_product(X, Y)
278 endif
279
280 nn = Dot_Product(n, n)
281
282 if (nn == 0) then
283 angle = 0.0_wp
284 else
285 n = n/sqrt(nn)
286 alpha = vector_product(n, X)
287
288 beta = X - Dot_Product(n, X) * n
289 gamma = Y - Dot_Product(n, Y) * n
290 angle = atan2(Dot_Product(alpha,gamma), Dot_Product(beta,gamma))
291 endif
292
293end function vector_angle
294
295
296function impact_parameter(r_leo, r_gns, bangle) result(impact)
297
298! 1.1 Declarations
299! ----------------
300
301 implicit none
302
303 real(wp), dimension(3), intent(in) :: r_leo ! LEO position vector (ECF)
304 real(wp), dimension(3), intent(in) :: r_gns ! GPS position vector (ECF)
305 real(wp), optional, intent(in) :: bangle ! Bending angle
306 real(wp) :: impact ! Impact parameter
307
308 real(wp) :: r0 ! Length of r_leo
309 real(wp) :: r1 ! Length of r_gns
310 real(wp) :: omega ! Complementary to r_leo^r_gns - bangle
311 real(wp) :: talpha ! Tan(r_leo^(r_leo-r_gns))
312
313! 1.2 Length of vectors r_leo and r_gns
314! -------------------------------------
315
316 r0 = Sqrt(Dot_Product(r_leo, r_leo))
317 r1 = Sqrt(Dot_Product(r_gns, r_gns))
318
319! 1.3 Find vector angle between r_leo and r_gns
320! ---------------------------------------------
321
322 omega = Pi - vector_angle(r_leo, r_gns)
323
324 if (present(bangle)) then
325 omega = omega + bangle
326 endif
327
328! 1.4 Determine impact parameter by trigonometry
329! ----------------------------------------------
330
331 talpha = r1*Sin(omega) / (r0 + r1*Cos(omega))
332
333 impact = r0*talpha / Sqrt(1.0_wp + talpha**2)
334
335end function impact_parameter
336
337
338
339
340end program azi
341
342
343
344
345
346
347
348
349