Ticket #457: ropp_pp_wopt_propagate_to_leo.f90

File ropp_pp_wopt_propagate_to_leo.f90, 22.9 KB (added by Ian Culverwell, 8 years ago)

ropp_pp_wopt_propagate_to_leo.f90

Line 
1! $Id: ropp_pp_wopt_propagate_to_leo.f90 r4887 2016-05-25 15:48:28Z sti $
2
3!****s* WaveOpticsPropagator/ropp_pp_wopt_propagate_to_leo *
4!
5! NAME
6! ropp_pp_wopt_propagate_to_leo
7!
8! SYNOPSIS
9!
10! CALL ropp_pp_wopt_propagate_to_leo(nsample, U0, k, ypos, Dy, x_plane, &
11! x_leo, y_leo, s_geom, U_leo, phase_leo, amp_leo)
12!
13! DESCRIPTION
14!
15! Propagate complex signal (phase and amplitude) from final screen to LEO
16!
17! INPUTS
18! INTEGER, INTENT(IN) :: nsample
19! COMPLEX(wp), INTENT(IN) :: U0
20! REAL(wp), INTENT(IN) :: k
21! REAL(wp), INTENT(IN) :: ypos
22! REAL(wp), INTENT(IN) :: Dy
23! REAL(wp), INTENT(IN) :: x_plane
24! REAL(wp), INTENT(IN) :: x_leo,y_leo
25! REAL(wp), INTENT(IN) :: s_geom
26!
27! OUTPUT
28! COMPLEX(wp), INTENT(OUT) :: U_leo
29! REAL(wp), INTENT(OUT) :: phase_leo, amp_leo
30!
31! NOTES
32! See RSR 27(?) for details
33!
34! AUTHOR
35! ECMWF, Reading, UK.
36! Any comments on this software should be given via the ROM SAF
37! Helpdesk at http://www.romsaf.org
38!
39! COPYRIGHT
40! (c) EUMETSAT. All rights reserved.
41! For further details please refer to the file COPYRIGHT
42! which you should have received as part of this distribution.
43!
44!****
45
46SUBROUTINE ropp_pp_wopt_propagate_to_leo(nsample, U0, k, ypos, Dy, x_plane, &
47 x_leo, y_leo, s_geom, U_leo, phase_leo, amp_leo)
48
49!-------------------------------------------------------------------------------
50! 1. Declarations
51!-------------------------------------------------------------------------------
52
53 USE typesizes, ONLY: wp => EightByteReal
54 USE ropp_pp_wopt, not_this => ropp_pp_wopt_propagate_to_leo
55 USE ropp_pp_constants, ONLY: pi1 => pi
56use ropp_io
57use ropp_io_types, only: roprof
58
59 IMPLICIT NONE
60
61! Input
62 INTEGER, INTENT(IN) :: nsample ! Number of points used in fitting
63 COMPLEX(wp), INTENT(IN) :: U0(:) ! Complex field in source plane
64 ! (SIZE(U0) must be a power of 2)
65 REAL(wp), INTENT(IN) :: k ! Wave number (2*Pi/Wavelength)
66 REAL(wp), INTENT(IN) :: Dy ! Discretization step in source plane
67 REAL(wp), INTENT(IN) :: ypos(:) ! y of screen
68 REAL(wp), INTENT(IN) :: x_plane ! x coordinate of source plane
69 REAL(wp), INTENT(IN) :: x_leo(:), y_leo(:) ! x,y of LEO orbit
70 REAL(wp), INTENT(IN) :: s_geom(:) ! GPS to LEO distance
71
72! Output
73 COMPLEX(wp), INTENT(OUT) :: U_leo(SIZE(x_leo)) ! Propagated field at LEO
74 REAL(wp), INTENT(OUT) :: phase_leo(SIZE(x_leo)), amp_leo(SIZE(x_leo)) ! phase and amp at LEO
75
76! Local
77 COMPLEX(wp), PARAMETER :: Ci = (0.0_wp, 1.0_wp) ! Sqrt(-1)
78
79 INTEGER :: ny ! Number of points on each screen
80 INTEGER :: n_leo ! Number of points in LEO orbit
81 INTEGER :: ny_av, ny_av2 ! Number of data points
82 INTEGER :: i ! Array index
83 INTEGER :: j ! Index for LEO position
84 INTEGER :: kk ! Loop index
85 INTEGER :: imin, imax ! Sub array bounds
86
87 REAL(wp), PARAMETER :: amplitude_cutoff = 1.0E-6_wp ! amplitudes < amplitude_cutoff*max_amp are ignored
88 REAL(wp) :: phase(SIZE(ypos)), amplitude(SIZE(ypos))
89 REAL(wp) :: amp_sum
90 REAL(wp) :: Dx
91 REAL(wp) :: twopi
92
93 COMPLEX(wp) :: U(SIZE(ypos)) ! For storing the arrays
94
95 COMPLEX(wp), ALLOCATABLE :: U_av(:) ! For computing integral
96
97 REAL(wp), ALLOCATABLE :: y_av(:), phase_grad(:), phase_start(:)
98 REAL(wp), ALLOCATABLE :: amp_grad(:), amp_start(:)
99
100 REAL(wp) :: c_fresnel(2), s_fresnel(2)
101 REAL(wp) :: Aconst, ymin, ymax, Amp_max
102
103 REAL(wp) :: aval,bval,cval
104
105 REAL(wp) :: y_upper,y_lower
106 REAL(wp) :: U_upper,U_lower
107 REAL(wp) :: real_term1,imag_term1
108 REAL(wp) :: real_term2,imag_term2
109 REAL(wp) :: rnsample
110
111 real(wp), allocatable :: amp_recon(:), amp_recon2(:)
112 real(wp), allocatable :: phase_recon(:), phase_recon2(:)
113 real(wp), allocatable :: xdata(:, :), ydata(:, :), ydata2(:, :), ydataout(:, :), bigxdata(:), bigydata(:)
114 REAL(wp), ALLOCATABLE :: amp_grad2(:), amp_start2(:)
115 REAL(wp), ALLOCATABLE :: phase_grad2(:), phase_start2(:)
116
117 type(roprof) :: diag
118!-------------------------------------------------------------------------------
119! 2. Initialise variables
120!-------------------------------------------------------------------------------
121
122 twopi = 2.0_wp * pi1
123
124 ny = SIZE(ypos)
125print*,'ny = ', ny
126 n_leo = SIZE(x_leo)
127
128 rnsample = REAL(nsample, KIND=wp)
129print*,'nsample = ', nsample
130
131 U(:) = U0(:) ! complex amplitude at final screen
132
133!-------------------------------------------------------------------------------
134! 3. Split the final screen into a set of vertical intervals
135!-------------------------------------------------------------------------------
136
137! 3.1 estimate the geometry
138! -------------------------
139
140 ny_av = INT( ny/nsample )
141print*,'ny_av = ', ny_av
142
143! 3.2 find the maximum amplitude at final screen
144! ----------------------------------------------
145
146 Amp_max = MAXVAL(ABS(U(:)))
147
148! 3.3 truncate summation on final screen
149! --------------------------------------
150
151 ymin = MINVAL( ypos, MASK = ABS(U) > amplitude_cutoff*Amp_max )
152
153 ymax = MAXVAL( ypos, MASK = ABS(U) > amplitude_cutoff*Amp_max )
154print*,'ymin = ', ymin
155print*,'ymax = ', ymax
156print*,'ypos(1) = ', ypos(1)
157print*,'ypos(ny) = ', ypos(ny)
158
159!-------------------------------------------------------------------------------
160! 4. Decompose complex signal into amplitude and accumulated phase
161!-------------------------------------------------------------------------------
162
163 CALL ropp_pp_wopt_phase_and_amplitude(U, phase, amplitude)
164
165!-------------------------------------------------------------------------------
166! 5. Allocate the arrays
167!-------------------------------------------------------------------------------
168
169 ny_av2 = 0
170
171 DO i = 1,ny_av
172
173 IF (ypos((i-1)*nsample+1) > ymax) EXIT
174
175 IF (ypos((i-1)*nsample+1) > ymin) ny_av2 = ny_av2 + 1
176
177 ENDDO
178
179print*,'ny_av2 = ', ny_av2
180
181 ALLOCATE (U_av(ny_av2))
182 ALLOCATE (phase_grad(ny_av2))
183 ALLOCATE (phase_start(ny_av2))
184 ALLOCATE (amp_grad(ny_av2))
185 ALLOCATE (amp_start(ny_av2))
186 ALLOCATE (y_av(ny_av2+1))
187 ALLOCATE (amp_recon(ny)) ; amp_recon = -1.0_wp
188 ALLOCATE (phase_recon(ny)) ; phase_recon = -1.0_wp
189 ALLOCATE (amp_recon2(ny)) ; amp_recon2 = -1.0_wp
190 ALLOCATE (phase_recon2(ny)) ; phase_recon2 = -1.0_wp
191 ALLOCATE (amp_grad2(ny_av2), amp_start2(ny_av2))
192 ALLOCATE (phase_grad2(ny_av2), phase_start2(ny_av2))
193
194! initialise
195
196 U_av(:) = (0.0_wp, 0.0_wp)
197
198 kk = 1
199
200!-------------------------------------------------------------------------------
201! 6. Make linear fit to phase and amplitude over each vertical interval
202!-------------------------------------------------------------------------------
203
204 DO i = 1,ny_av
205
206 IF (kk == ny_av2+1) EXIT
207
208 imin = (i-1)*nsample + 1 ; imax = imin + nsample - 1
209
210 IF (ypos(imin) > ymin) THEN
211
212 y_av(kk) = ypos(imin)
213
214 CALL lsq_fit( (/ (REAL(j, KIND=wp)*dy, j=0,nsample-1) /), &
215 phase(imin:imax), phase_grad(kk), phase_start(kk) )
216
217 CALL lsq_fit( (/ (REAL(j, KIND=wp)*dy, j=0,nsample-1) /), &
218 amplitude(imin:imax), amp_grad(kk), amp_start(kk) )
219
220 amp_recon(imin:imax) = amp_start(kk) + &
221 amp_grad(kk) * (/ (REAL(j, KIND=wp)*dy, j=0,nsample-1) /)
222
223 phase_recon(imin:imax) = phase_start(kk) + &
224 phase_grad(kk) * (/ (REAL(j, KIND=wp)*dy, j=0,nsample-1) /)
225
226 U_av(kk) = CMPLX(COS(phase_start(kk)), SIN(phase_start(kk)), KIND=wp) ! phase of signal at kk
227
228 kk = kk + 1 ! next vertical interval
229
230 ENDIF
231
232 ENDDO
233
234print*,'Overall (amp-amp_fit) sum of squares = ', SUM((amplitude-amp_recon)**2, MASK=(ypos >= ymin) .AND. (ypos <= ymax))
235print*,'Overall (phase-phase_fit) sum of squares = ', SUM((phase-phase_recon)**2, MASK=(ypos >= ymin) .AND. (ypos <= ymax))
236
237! Try continuous lsq fit
238
239 ALLOCATE ( xdata(nsample, ny_av2), ydata(nsample, ny_av2), ydata2(nsample, ny_av2), ydataout(nsample, ny_av2), bigxdata(ny_av2+1), bigydata(ny_av2+1) )
240 xdata = -1.0_wp ; ydata = -1.0_wp ; ydata2 = -1.0_wp ; ydataout = -1.0_wp ; bigxdata = -1.0_wp ; bigydata = -1.0_wp
241
242 kk = 1
243 DO i = 1,ny_av
244
245 IF (kk == ny_av2+1) EXIT
246
247 imin = (i-1)*nsample + 1 ; imax = imin + nsample - 1
248
249 IF (ypos(imin) > ymin) THEN
250
251 xdata(:, kk) = ypos(imin:imax)
252 ydata(:, kk) = amplitude(imin:imax)
253 ydata2(:, kk) = phase(imin:imax)
254
255 bigxdata(kk) = ypos(imin)
256
257 kk = kk + 1 ! next vertical interval
258
259 ENDIF
260
261 ENDDO
262
263 bigxdata(kk) = ypos(imax+1) ! kk = ny_av2 + 1
264 bigydata = bigxdata
265
266! Amps
267 CALL cont_lsq_fit (xdata, ydata, bigxdata, bigydata, ydataout)
268
269 kk = 1
270 DO i = 1,ny_av
271
272 IF (kk == ny_av2+1) EXIT
273
274 imin = (i-1)*nsample + 1 ; imax = imin + nsample - 1
275
276 IF (ypos(imin) > ymin) THEN
277 amp_grad2(kk) = (bigydata(kk+1) - bigydata(kk)) / (bigxdata(kk+1) - bigxdata(kk))
278 amp_start2(kk) = (bigxdata(kk+1)*bigydata(kk) - bigxdata(kk)*bigydata(kk+1)) / (bigxdata(kk+1) - bigxdata(kk))
279 amp_start2(kk) = amp_start2(kk) + amp_grad2(kk)*ypos(imin) ! Need to subtract offset
280 amp_recon2(imin:imax) = ydataout(:, kk)
281!!!!!!!!!!!!
282 amp_grad(kk) = amp_grad2(kk) ; amp_start(kk) = amp_start2(kk)
283!!!!!!!!!!!!
284 kk = kk + 1 ! next vertical interval
285 END IF
286
287 ENDDO
288
289print*,'Overall (amp-amp_fit2) sum of squares = ', SUM((amplitude-amp_recon2)**2, MASK=(ypos >= ymin) .AND. (ypos <= ymax))
290print*,'Overall (amp_fit1-amp_fit2) sum of squares = ', SUM((amp_recon-amp_recon2)**2, MASK=(ypos >= ymin) .AND. (ypos <= ymax))
291
292! Phases
293 CALL cont_lsq_fit (xdata, ydata2, bigxdata, bigydata, ydataout)
294
295 kk = 1
296 DO i = 1,ny_av
297
298 IF (kk == ny_av2+1) EXIT
299
300 imin = (i-1)*nsample + 1 ; imax = imin + nsample - 1
301
302 IF (ypos(imin) > ymin) THEN
303 phase_grad2(kk) = (bigydata(kk+1) - bigydata(kk)) / (bigxdata(kk+1) - bigxdata(kk))
304 phase_start2(kk) = (bigxdata(kk+1)*bigydata(kk) - bigxdata(kk)*bigydata(kk+1)) / (bigxdata(kk+1) - bigxdata(kk))
305 phase_start2(kk) = phase_start2(kk) + phase_grad2(kk)*ypos(imin) ! Need to subtract offset
306 phase_recon2(imin:imax) = ydataout(:, kk)
307!!!!!!!!!!!!
308 phase_grad(kk) = phase_grad2(kk) ; phase_start(kk) = phase_start2(kk)
309!!!!!!!!!!!!
310 kk = kk + 1 ! next vertical interval
311 END IF
312
313 ENDDO
314
315print*,'Overall (phase-phase_fit2) sum of squares = ', SUM((phase-phase_recon2)**2, MASK=(ypos >= ymin) .AND. (ypos <= ymax))
316print*,'Overall (phase_fit1-phase_fit2) sum of squares = ', SUM((phase_recon-phase_recon2)**2, MASK=(ypos >= ymin) .AND. (ypos <= ymax))
317
318
319
320 DEALLOCATE ( bigydata, bigxdata, ydataout, ydata2, ydata, xdata )
321
322
323 CALL ropp_io_addvar_rodataD1d( diag, &
324 name = "Amp_grad", &
325 long_name = "Gradient of amplitude fits", &
326 units = "", &
327 range = (/ -1000000.0_wp, 1000000.0_wp /), &
328 DATA = amp_grad )
329
330 CALL ropp_io_addvar_rodataD1d( diag, &
331 name = "Amp_int", &
332 long_name = "Intercept of amplitude fits", &
333 units = "", &
334 range = (/ -1000000.0_wp, 1000000.0_wp /), &
335 DATA = amp_start )
336
337 CALL ropp_io_addvar_rodataD1d( diag, &
338 name = "Amp_grad2", &
339 long_name = "Gradient of amplitude fits2", &
340 units = "", &
341 range = (/ -1000000.0_wp, 1000000.0_wp /), &
342 DATA = amp_grad2 )
343
344 CALL ropp_io_addvar_rodataD1d( diag, &
345 name = "Amp_int2", &
346 long_name = "Intercept of amplitude fits2", &
347 units = "", &
348 range = (/ -1000000.0_wp, 1000000.0_wp /), &
349 DATA = amp_start2 )
350
351 WHERE ( (ypos < ymin) .or. (ypos > ymax) ) amplitude=-1.0_wp
352 CALL ropp_io_addvar_rodataD1d( diag, &
353 name = "Amp", &
354 long_name = "High res amplitude", &
355 units = "", &
356 range = (/ -1000000.0_wp, 1000000.0_wp /), &
357 DATA = amplitude )
358
359 WHERE ( (ypos < ymin) .or. (ypos > ymax) ) amp_recon=-1.0_wp
360 CALL ropp_io_addvar_rodataD1d( diag, &
361 name = "Amp_recon", &
362 long_name = "Reconstructed amplitude fits", &
363 units = "", &
364 range = (/ -1000000.0_wp, 1000000.0_wp /), &
365 DATA = amp_recon )
366
367 WHERE ( (ypos < ymin) .or. (ypos > ymax) ) amp_recon2=-1.0_wp
368 CALL ropp_io_addvar_rodataD1d( diag, &
369 name = "Amp_recon2", &
370 long_name = "Reconstructed amp fits (cts)", &
371 units = "", &
372 range = (/ -1000000.0_wp, 1000000.0_wp /), &
373 DATA = amp_recon2 )
374
375
376 CALL ropp_io_addvar_rodataD1d( diag, &
377 name = "Phase_grad", &
378 long_name = "Gradient of phase fits", &
379 units = "", &
380 range = (/ -1000000.0_wp, 1000000.0_wp /), &
381 DATA = phase_grad )
382
383 CALL ropp_io_addvar_rodataD1d( diag, &
384 name = "Phase_int", &
385 long_name = "Intercept of phase fits", &
386 units = "", &
387 range = (/ -1000000.0_wp, 1000000.0_wp /), &
388 DATA = phase_start )
389
390 CALL ropp_io_addvar_rodataD1d( diag, &
391 name = "Phase_grad2", &
392 long_name = "Gradient of phase fits2", &
393 units = "", &
394 range = (/ -1000000.0_wp, 1000000.0_wp /), &
395 DATA = phase_grad2 )
396
397 CALL ropp_io_addvar_rodataD1d( diag, &
398 name = "Phase_int2", &
399 long_name = "Intercept of phase fits2", &
400 units = "", &
401 range = (/ -1000000.0_wp, 1000000.0_wp /), &
402 DATA = phase_start2 )
403
404 WHERE ( (ypos < ymin) .or. (ypos > ymax) ) phase=-1.0_wp
405 CALL ropp_io_addvar_rodataD1d( diag, &
406 name = "Phase", &
407 long_name = "High res phase", &
408 units = "", &
409 range = (/ -1000000.0_wp, 1000000.0_wp /), &
410 DATA = phase )
411
412 WHERE ( (ypos < ymin) .or. (ypos > ymax) ) phase_recon=-1.0_wp
413 CALL ropp_io_addvar_rodataD1d( diag, &
414 name = "Phase_recon", &
415 long_name = "Reconstructed phase fits", &
416 units = "", &
417 range = (/ -1000000.0_wp, 1000000.0_wp /), &
418 DATA = phase_recon )
419
420 WHERE ( (ypos < ymin) .or. (ypos > ymax) ) phase_recon2=-1.0_wp
421 CALL ropp_io_addvar_rodataD1d( diag, &
422 name = "Phase_recon2", &
423 long_name = "Reconstructed phase fits (cts)", &
424 units = "", &
425 range = (/ -1000000.0_wp, 1000000.0_wp /), &
426 DATA = phase_recon2 )
427
428
429 CALL ropp_io_write(diag, 'data/diag.nc', append=.false., ranchk=.false.)
430
431 CALL ropp_io_free(diag)
432
433 DEALLOCATE (amp_recon, phase_recon, amp_recon2, phase_recon2)
434
435! uppermost height
436
437 y_av(ny_av2+1) = y_av(ny_av2) + 1.0E20_wp ! force upper limit to ~infinity
438
439!-------------------------------------------------------------------------------
440! 7. Compute the complex amplitude at each LEO position
441!-------------------------------------------------------------------------------
442
443! 7.1 initialise the field at LEO
444! -------------------------------
445
446 U_leo(:) = (0.0_wp, 0.0_wp)
447
448! 7.2 loop through the observations
449! ---------------------------------
450
451 DO i = 1,n_leo
452
453 Dx = x_leo(i) - x_plane ! distance from plane to LEO at position
454
455! for computing beta
456
457 aval = 0.5_wp*k/Dx
458
459 Aconst = SQRT(2.0_wp/(aval*pi1))
460
461! 7.3 loop through the vertical integration intervals
462! ---------------------------------------------------
463
464 DO j = 1, ny_av2
465
466 bval = 0.5_wp*phase_grad(j)
467
468 cval = phase_grad(j)*(y_leo(i)-y_av(j))
469
470! 7.3.1 the part coming from the constant term in the amplitude
471! -------------------------------------------------------------
472
473! limits of integration
474
475 y_upper = Aconst*(aval*(y_av(j+1)-y_leo(i))+bval)
476
477 y_lower = Aconst*(aval*(y_av(j )-y_leo(i))+bval)
478
479 IF (j == 1) y_lower = y_lower - 1.0E20_wp ! force the lower limit to ~infinity
480
481 CALL ropp_pp_wopt_fresnel((/y_upper, y_lower/), s_fresnel, c_fresnel)
482
483 amp_sum = amp_start(j) + amp_grad(j)*(y_leo(i) - y_av(j) - bval/aval)
484
485 real_term1 = SQRT(0.5_wp*Pi1/aval)*amp_sum*(c_fresnel(1)-c_fresnel(2))
486
487 imag_term1 = SQRT(0.5_wp*Pi1/aval)*amp_sum*(s_fresnel(1)-s_fresnel(2))
488
489! 7.3.2 the part coming from the linearly varying term in the amplitude
490! ---------------------------------------------------------------------
491
492 real_term2 = 0.0_wp
493 imag_term2 = 0.0_wp
494
495 IF ( j > 1 .AND. j < ny_av2 ) THEN ! Assume constant extrapolation in first and last samples (~ amp_grad = 0)
496
497! limits of integration
498 U_upper = (aval*(y_av(j+1)-y_leo(i))+bval)**2/aval
499
500 U_lower = (aval*(y_av(j )-y_leo(i))+bval)**2/aval
501
502 real_term2 = (0.5_wp*amp_grad(j)/aval)*(SIN(U_upper)-SIN(U_lower))
503
504 imag_term2 = -(0.5_wp*amp_grad(j)/aval)*(COS(U_upper)-COS(U_lower))
505
506 END IF
507
508! 7.3.3 combine them
509! ------------------
510
511 U_leo(i) = U_leo(i) + U_av(j) * &
512 CMPLX(real_term1+real_term2, imag_term1+imag_term2, KIND=wp) * &
513 EXP(Ci*MODULO((aval*cval-bval**2)/aval, twopi))
514
515 ENDDO
516
517! 7.4 add k.Dx factor and subtract the vacuum delay, s_geom
518! -----------------------------------------------------------
519
520 U_leo(i) = U_leo(i) * &
521 EXP(Ci*MODULO(k*(Dx-s_geom(i)), twopi)) * &
522 SQRT(k/(twopi*Dx))
523
524 ENDDO
525
526!-------------------------------------------------------------------------------
527! 8. Compute the amplitude and accumulated phase at LEO
528!-------------------------------------------------------------------------------
529
530 CALL ropp_pp_wopt_phase_and_amplitude_LEO(U_leo, phase_leo, amp_leo)
531
532! convert to metres by dividing by wave number
533
534 phase_leo(:) = phase_leo(:) / k
535
536!-------------------------------------------------------------------------------
537! 9. Tidy up
538!-------------------------------------------------------------------------------
539
540 DEALLOCATE(y_av)
541 DEALLOCATE(amp_start)
542 DEALLOCATE(amp_grad)
543 DEALLOCATE(phase_start)
544 DEALLOCATE(phase_grad)
545 DEALLOCATE(U_av)
546
547 DEALLOCATE (amp_grad2, amp_start2)
548
549
550CONTAINS
551
552!-------------------------------------------------------------------------------
553! 10. Standard least-squares fit to y = m*x + c
554!-------------------------------------------------------------------------------
555
556 SUBROUTINE lsq_fit (x, y, m, c)
557
558 REAL(wp), DIMENSION(:), INTENT(IN) :: x, y
559 REAL(wp), INTENT(OUT) :: m, c
560 REAL(wp) :: rnx
561 REAL(wp) :: xbar, ybar, xxbar, xybar
562
563 rnx = REAL(SIZE(x), KIND=wp)
564
565 xbar = SUM(x) / rnx
566
567 ybar = SUM(y) / rnx
568
569 xxbar = SUM(x*x) / rnx
570
571 xybar = SUM(x*y) / rnx
572
573 m = (xybar - xbar*ybar) / (xxbar - xbar*xbar)
574
575 c = ybar - m * xbar
576
577 END SUBROUTINE lsq_fit
578
579!-------------------------------------------------------------------------------
580! 11. Continuous piecewise-linear least-squares fit
581!-------------------------------------------------------------------------------
582
583 SUBROUTINE cont_lsq_fit (x, y, bigx, bigy, yout)
584
585! Make piecewise linear LSQ fit to {{(x^i_J, y^i_J), i=1 ... n_J}, J=1, N)
586! by constructing the ordinates Y_J at the given abscissae X_J.
587! x^i_J is held in a 2D array x(i, J).
588! Currently we assume n_J is constant.
589
590! I/O
591 REAL(wp), INTENT(IN) :: x(:, :), y(:, :), bigx(:)
592 REAL(wp), INTENT(INOUT) :: bigy(:)
593 REAL(wp), OPTIONAL, INTENT(INOUT) :: yout(:, :)
594! Local
595 INTEGER :: nj, bign, bigj
596 REAL(wp), ALLOCATABLE :: xbar(:), ybar(:), xxbar(:), xybar(:)
597 REAL(wp), ALLOCATABLE :: dl(:), dd(:), du(:), rhs(:)
598 REAL(wp) :: m, c
599
600 nj = SIZE(x, 1)
601 bign = SIZE(x, 2)
602
603! Calculate means
604
605 ALLOCATE ( xbar(1:bign), xxbar(1:bign), ybar(1:bign), xybar(1:bign) )
606
607 xbar = SUM(x, DIM=1) / nj
608 xxbar = SUM(x*x, DIM=1) / nj
609 ybar = SUM(y, DIM=1) / nj
610 xybar = SUM(x*y, DIM=1) / nj
611
612! Define tridiagonal matrix
613
614 ALLOCATE ( dl(1:bign+1), dd(1:bign+1), du(1:bign+1), rhs(1:bign+1) )
615
616! Top row of tridiag
617 bigj = 1
618 dl(bigj) = 0.0_wp
619 dd(bigj) = xxbar(bigj) - xbar(bigj)*2.0_wp*bigx(bigj+1) + bigx(bigj+1)**2
620 du(bigj) = -(xxbar(bigj) - xbar(bigj)*(bigx(bigj)+bigx(bigj+1)) + bigx(bigj)*bigx(bigj+1))
621 rhs(bigj) = -(bigx(bigj+1)-bigx(bigj)) * (xybar(bigj) - ybar(bigj)*bigx(bigj+1))
622
623! Bottom row of tridiag
624 bigj = bign + 1
625 dl(bigj) = -(xxbar(bigj-1) - xbar(bigj-1)*(bigx(bigj)+bigx(bigj-1)) + bigx(bigj)*bigx(bigj-1))
626 dd(bigj) = xxbar(bigj-1) - xbar(bigj-1)*2.0_wp*bigx(bigj-1) + bigx(bigj-1)**2
627 du(bigj) = 0.0_wp
628 rhs(bigj) = (bigx(bigj)-bigx(bigj-1)) * (xybar(bigj-1) - ybar(bigj-1)*bigx(bigj-1))
629
630! Other rows of tridiag
631 DO bigj=2,bign
632 dl(bigj) = -(xxbar(bigj-1) - xbar(bigj-1)*(bigx(bigj)+bigx(bigj-1)) + bigx(bigj)*bigx(bigj-1))
633 dd(bigj) = xxbar(bigj-1) - xbar(bigj-1)*2.0_wp*bigx(bigj-1) + bigx(bigj-1)**2 + &
634 xxbar(bigj) - xbar(bigj)*2.0_wp*bigx(bigj+1) + bigx(bigj+1)**2
635 du(bigj) = -(xxbar(bigj) - xbar(bigj)*(bigx(bigj)+bigx(bigj+1)) + bigx(bigj)*bigx(bigj+1))
636 rhs(bigj) = (bigx(bigj)-bigx(bigj-1)) * (xybar(bigj-1) - ybar(bigj-1)*bigx(bigj-1)) - &
637 (bigx(bigj+1)-bigx(bigj)) * (xybar(bigj) - ybar(bigj)*bigx(bigj+1))
638 END DO
639
640! Tridiagonal inversion
641
642 CALL TRISOL(dl, dd, du, rhs, bigy)
643
644! Generate fit at original x values
645
646 IF ( PRESENT(yout) ) THEN
647
648 DO bigj=1,bign
649
650 m = (bigy(bigj+1) - bigy(bigj)) / (bigx(bigj+1) - bigx(bigj))
651 c = bigy(bigj+1) - m * bigx(bigj+1)
652 yout(:, bigj) = m * x(:, bigj) + c
653
654 END DO
655
656 END IF
657
658! Clean up
659
660 DEALLOCATE ( rhs, du, dd, dl )
661 DEALLOCATE (xybar, ybar, xxbar, xbar )
662
663 END SUBROUTINE cont_lsq_fit
664
665!-------------------------------------------------------------------------------
666! 12. Tridiagonal inversion
667!-------------------------------------------------------------------------------
668
669 SUBROUTINE TRISOL(a, b, c, d, x)
670
671! Solves
672!
673! ( b1 c1 ) (x1 ) = (d1 )
674! ( a2 b2 c2 ) (x2 ) = (d2 )
675! ( ai bi ci ) (xi ) = (di )
676! ( an-1 bn-1 cn-1) (xn-1) = (dn-1)
677! ( an bn) (xn ) = (dn )
678!
679! by Gaussian elimination without pivoting, eg
680! https://en.wikipedia.org/w/index.php?title=Tridiagonal_matrix_algorithm.
681! a1 and cn must be supplied (i.e. all vectors must have dimension (1:n)),
682! even though they are not used.
683
684! I/O
685 REAL(wp), INTENT(IN) :: a(:), b(:), c(:), d(:)
686 REAL(wp), INTENT(INOUT) :: x(:)
687! Local
688 INTEGER :: ii, nn
689 REAL(wp), ALLOCATABLE :: c1(:), d1(:)
690
691 nn = SIZE(a)
692 ALLOCATE (c1(nn), d1(nn))
693
694 c1(1) = c(1) / b(1)
695 DO ii=2,nn-1
696 c1(ii) = c(ii) / (b(ii) - a(ii)*c1(ii-1))
697 END DO
698
699 d1(1) = d(1) / b(1)
700 DO ii=2,nn
701 d1(ii) = (d(ii) - a(ii)*d1(ii-1)) / (b(ii) - a(ii)*c1(ii-1))
702 END DO
703
704 x(nn) = d1(nn)
705 DO ii=nn-1,1,-1
706 x(ii) = d1(ii) - c1(ii)*x(ii+1)
707 END DO
708
709 DEALLOCATE (d1, c1)
710
711 END SUBROUTINE TRISOL
712
713END SUBROUTINE ropp_pp_wopt_propagate_to_leo