! $Id: ropp_pp_wopt_propagate_to_leo.f90 r4887 2016-05-25 15:48:28Z sti $ !****s* WaveOpticsPropagator/ropp_pp_wopt_propagate_to_leo * ! ! NAME ! ropp_pp_wopt_propagate_to_leo ! ! SYNOPSIS ! ! CALL ropp_pp_wopt_propagate_to_leo(nsample, U0, k, ypos, Dy, x_plane, & ! x_leo, y_leo, s_geom, U_leo, phase_leo, amp_leo) ! ! DESCRIPTION ! ! Propagate complex signal (phase and amplitude) from final screen to LEO ! ! INPUTS ! INTEGER, INTENT(IN) :: nsample ! COMPLEX(wp), INTENT(IN) :: U0 ! REAL(wp), INTENT(IN) :: k ! REAL(wp), INTENT(IN) :: ypos ! REAL(wp), INTENT(IN) :: Dy ! REAL(wp), INTENT(IN) :: x_plane ! REAL(wp), INTENT(IN) :: x_leo,y_leo ! REAL(wp), INTENT(IN) :: s_geom ! ! OUTPUT ! COMPLEX(wp), INTENT(OUT) :: U_leo ! REAL(wp), INTENT(OUT) :: phase_leo, amp_leo ! ! REFERENCES ! Description of wave optics modelling in ROPP-9 and suggested ! improvements for ROPP-9.1, S. B. Healy, ROM SAF Report 28, ! SAF/ROM/ECMWF/REP/RSR/028 ! ! AUTHOR ! ECMWF, Reading, UK. ! Any comments on this software should be given via the ROM SAF ! Helpdesk at http://www.romsaf.org ! ! COPYRIGHT ! (c) EUMETSAT. All rights reserved. ! For further details please refer to the file COPYRIGHT ! which you should have received as part of this distribution. ! !**** SUBROUTINE ropp_pp_wopt_propagate_to_leo(nsample, U0, k, ypos, Dy, x_plane, & x_leo, y_leo, s_geom, U_leo, phase_leo, amp_leo) !------------------------------------------------------------------------------- ! 1. Declarations !------------------------------------------------------------------------------- USE typesizes, ONLY: wp => EightByteReal USE ropp_pp_wopt, not_this => ropp_pp_wopt_propagate_to_leo USE ropp_pp_constants, ONLY: pi1 => pi IMPLICIT NONE ! Input INTEGER, INTENT(IN) :: nsample ! number of points used in fitting COMPLEX(wp), INTENT(IN) :: U0(:) ! Complex field in source plane ! (SIZE(U0) must be a power of 2) REAL(wp), INTENT(IN) :: k ! Wave number (2*Pi/Wavelength) REAL(wp), INTENT(IN) :: Dy ! Discretization step in source plane REAL(wp), INTENT(IN) :: ypos(:) ! y of screen REAL(wp), INTENT(IN) :: x_plane ! x coordinate of in source plane REAL(wp), INTENT(IN) :: x_leo(:), y_leo(:) ! x,y of LEO orbit REAL(wp), INTENT(IN) :: s_geom(:) ! GPS to LEO distance ! Output COMPLEX(wp), INTENT(OUT) :: U_leo(SIZE(x_leo)) ! Propagated field at LEO REAL(wp), INTENT(OUT) :: phase_leo(SIZE(x_leo)), amp_leo(SIZE(x_leo)) ! phase and amp at LEO ! Local COMPLEX(wp), PARAMETER :: Ci = (0.0_wp, 1.0_wp) ! Sqrt(-1) INTEGER :: ny ! number of points on each screen INTEGER :: n_leo ! number of points in LEO orbit INTEGER :: ny_av, ny_av2 ! Number of data points INTEGER :: i ! Array index INTEGER :: j ! index for LEO position INTEGER :: kk ! loop index INTEGER :: imin, imax ! sub array bounds REAL(wp), PARAMETER :: amplitude_cutoff = 1.0E-6_wp ! amplitudes < amplitude_cutoff*max_amp are ignored REAL(wp) :: phase(SIZE(ypos)), amplitude(SIZE(ypos)) REAL(wp) :: amp_sum REAL(wp) :: Dx REAL(wp) :: twopi COMPLEX(wp) :: U(SIZE(ypos)) ! for storing the arrays COMPLEX(wp), ALLOCATABLE :: U_av(:) ! for computing integral REAL(wp), ALLOCATABLE :: y_av(:), y_sample(:), phase_grad(:), phase_start(:) REAL(wp), ALLOCATABLE :: amp_grad(:), amp_start(:) REAL(wp) :: c_fresnel(2), s_fresnel(2) REAL(wp) :: Aconst, ymin, ymax, Amp_max REAL(wp) :: aval,bval,cval REAL(wp) :: y_upper,y_lower REAL(wp) :: U_upper,U_lower REAL(wp) :: real_term1,imag_term1 REAL(wp) :: real_term2,imag_term2 REAL(wp) :: ymid,rad_residual !------------------------------------------------------------------------------- ! 2. Initialise variables !------------------------------------------------------------------------------- twopi = 2.0_wp * pi1 ny = SIZE(ypos) n_leo = SIZE(x_leo) U(:) = U0(:) ! complex amplitude at final screen !------------------------------------------------------------------------------- ! 3. Split the final screen into a set of vertical intervals !------------------------------------------------------------------------------- ! 3.1 estimate the geometry ! ------------------------- ny_av = INT( ny / nsample ) ! 3.2 find the maximum amplitude at final screen ! ---------------------------------------------- Amp_max = MAXVAL(ABS(U(:))) ! 3.3 truncate summation on final screen ! -------------------------------------- ymin = MINVAL( ypos, MASK = ABS(U) > amplitude_cutoff*Amp_max ) ymax = MAXVAL( ypos, MASK = ABS(U) > amplitude_cutoff*Amp_max ) !------------------------------------------------------------------------------- ! 4. Decompose complex signal into amplitude and accumulated phase !------------------------------------------------------------------------------- CALL ropp_pp_wopt_phase_and_amplitude(U, phase, amplitude) !------------------------------------------------------------------------------- ! 5. Allocate the arrays !------------------------------------------------------------------------------- ny_av2 = 0 DO i = 1,ny_av IF (ypos((i-1)*nsample+1) > ymax) EXIT IF (ypos((i-1)*nsample+1) > ymin) ny_av2 = ny_av2 + 1 ENDDO ALLOCATE (U_av(ny_av2)) ALLOCATE (phase_grad(ny_av2)) ALLOCATE (phase_start(ny_av2)) ALLOCATE (amp_grad(ny_av2)) ALLOCATE (amp_start(ny_av2)) ALLOCATE (y_av(ny_av2+1)) ALLOCATE (y_sample(nsample)) ! initialise U_av(:) = (0.0_wp, 0.0_wp) y_sample = (/ (REAL(j, KIND=wp), j=0,nsample-1) /) * dy ! dummy coordinate across each screen kk = 1 !------------------------------------------------------------------------------- ! 6. Make linear fit to phase and amplitude over each vertical interval !------------------------------------------------------------------------------- DO i = 1,ny_av IF (kk == ny_av2+1) EXIT imin = (i-1)*nsample + 1 ; imax = imin + nsample - 1 IF (ypos(imin) > ymin) THEN y_av(kk) = ypos(imin) CALL lsq_fit( y_sample, phase(imin:imax), phase_grad(kk), phase_start(kk) ) CALL lsq_fit( y_sample, amplitude(imin:imax), amp_grad(kk), amp_start(kk) ) U_av(kk) = CMPLX(COS(phase_start(kk)), SIN(phase_start(kk)), KIND=wp) ! phase of signal at kk kk = kk + 1 ! next vertical interval ENDIF ENDDO ! uppermost height y_av(ny_av2+1) = y_av(ny_av2) + 1.0E20_wp ! force upper limit to ~infinity !------------------------------------------------------------------------------- ! 7. Compute the complex amplitude at each LEO position !------------------------------------------------------------------------------- ! 7.1 initialise the field at LEO ! ------------------------------- U_leo(:) = (0.0_wp, 0.0_wp) ! 7.2 loop through the observations ! --------------------------------- DO i = 1,n_leo Dx = x_leo(i) - x_plane ! distance from plane to LEO at position ! for computing beta aval = 0.5_wp*k/Dx Aconst = SQRT(2.0_wp/(aval*pi1)) ! 7.3 loop through the vertical integration intervals ! --------------------------------------------------- DO j = 1, ny_av2 bval = 0.5_wp*phase_grad(j) cval = phase_grad(j)*(y_leo(i)-y_av(j)) ! 7.3.1 the part coming from the constant term in the amplitude ! ------------------------------------------------------------- ! limits of integration y_upper = Aconst*(aval*(y_av(j+1)-y_leo(i))+bval) y_lower = Aconst*(aval*(y_av(j )-y_leo(i))+bval) IF (j == 1) y_lower = y_lower - 1.0E20_wp ! force the lower limit to ~infinity CALL ropp_pp_wopt_fresnel((/y_upper, y_lower/), s_fresnel, c_fresnel) amp_sum = amp_start(j) + amp_grad(j)*(y_leo(i) - y_av(j) - bval/aval) real_term1 = SQRT(0.5_wp*Pi1/aval)*amp_sum*(c_fresnel(1)-c_fresnel(2)) imag_term1 = SQRT(0.5_wp*Pi1/aval)*amp_sum*(s_fresnel(1)-s_fresnel(2)) ! 7.3.2 the part coming from the linearly varying term in the amplitude ! --------------------------------------------------------------------- real_term2 = 0.0_wp imag_term2 = 0.0_wp IF ( j > 1 .AND. j < ny_av2 ) THEN ! Assume constant extrapolation in first and last samples (~ amp_grad = 0) ! limits of integration U_upper = (aval*(y_av(j+1)-y_leo(i))+bval)**2/aval U_lower = (aval*(y_av(j )-y_leo(i))+bval)**2/aval real_term2 = (0.5_wp*amp_grad(j)/aval)*(SIN(U_upper)-SIN(U_lower)) imag_term2 = -(0.5_wp*amp_grad(j)/aval)*(COS(U_upper)-COS(U_lower)) END IF ! Approximate next term in the Fresnel expansion. Only include variation between screens. IF (j == 1 .OR. j == ny_av2) THEN ymid = y_av(j) ELSE ymid = 0.5_wp * (y_av(j) + y_av(j+1)) END IF rad_residual = -0.125_wp*k*(ymid - y_leo(i))**4 / Dx**3 ! 7.3.3 combine them ! ------------------ U_leo(i) = U_leo(i) + & U_av(j) * & CMPLX(real_term1+real_term2, imag_term1+imag_term2, KIND=wp) * & EXP(Ci*MODULO(((aval*cval-bval**2)/aval) + rad_residual, twopi)) END DO ! 7.4 add k.Dx factor and subtract the vacuum delay, s_geom ! ----------------------------------------------------------- U_leo(i) = U_leo(i) * SQRT(k/twopi/Dx) * & EXP(Ci*MODULO(k*(Dx-s_geom(i)), twopi)) END DO !------------------------------------------------------------------------------- ! 8. Compute the amplitude and accumulated phase at LEO !------------------------------------------------------------------------------- CALL ropp_pp_wopt_phase_and_amplitude_LEO(U_leo, phase_leo, amp_leo) ! convert to metres by dividing by wave number phase_leo(:) = phase_leo(:) / k !------------------------------------------------------------------------------- ! 9. Tidy up !------------------------------------------------------------------------------- DEALLOCATE(y_sample) DEALLOCATE(y_av) DEALLOCATE(amp_start) DEALLOCATE(amp_grad) DEALLOCATE(phase_start) DEALLOCATE(phase_grad) DEALLOCATE(U_av) CONTAINS !------------------------------------------------------------------------------- ! 10. Standard least-squares fit to y = m*x + c !------------------------------------------------------------------------------- SUBROUTINE lsq_fit (x, y, m, c) REAL(wp), DIMENSION(:), INTENT(IN) :: x, y REAL(wp), INTENT(OUT) :: m, c REAL(wp) :: rnx REAL(wp) :: xbar, ybar, xxbar, xybar rnx = REAL(SIZE(x), KIND=wp) xbar = SUM(x) / rnx ybar = SUM(y) / rnx xxbar = SUM(x*x) / rnx xybar = SUM(x*y) / rnx m = (xybar - xbar*ybar) / (xxbar - xbar*xbar) c = ybar - m * xbar END SUBROUTINE lsq_fit END SUBROUTINE ropp_pp_wopt_propagate_to_leo