! $Id: fsi_leo.f90 2022-11-02 09:51:28Z ychen $ !****s* bangle/fsi_leo * ! ! NAME ! fsi_leo - calculate the bending angle and impact height using FSI method ! ! SYNOPSIS ! call fsi_leo(ang, full_phase, amp, r1, r2, pp, ee, aa, pmin, pmax, apply_SH, nout) ! ! DESCRIPTION ! This subroutine calculate the bending angle and impact height using FSI method ! ! INPUTS ! REAL(wp), DIMENSION(:) :: ang central angle between 1st and 2nd radius-vectors ! REAL(wp), DIMENSION(:,:) :: full_phase full phase (m) after propagating to circular orbit[ch, t] ! REAL(wp), DIMENSION(:,:) :: amp Amplitudes [ch, time] ! REAL(wp) :: r1 Radius of GNSS (m) ! REAL(wp) :: r2 Radius of LEO (m) ! REAL(wp) :: pmin min impact height ! REAL(wp) :: pmax max impact height ! LOGICAL :: apply_SH flag to apply shadow height determination ! ! OUTPUT ! REAL(wp), DIMENSION(:,:) :: pp Impact parameters (m) ! REAL(wp), DIMENSION(:,:) :: ee Bending angles (rad) ! REAL(wp), DIMENSION(:,:) :: aa FSI amplitude ! INTEGER, DIMENSION(:), :: nout array size for impact and bending ! ! AUTHOR ! Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov ! Loknath Adhikari, UMD/CISESS, loknath.adhikari@noaa.gov ! ! COPYRIGHT ! Copyright (c) 2022-2023 Yong Chen ! For further details please refer to the file COPYRIGHT ! which you should have received as part of this distribution. ! !**** SUBROUTINE fsi_leo(ang, full_phase, amp, r1, r2, pp, ee, aa, pmin, pmax, apply_SH, nout) USE typesizes, ONLY: wp => EightByteReal USE ropp_pp_constants, ONLY: c_light, pi, f_L1, f_L2 USE ropp_pp, ONLY: ropp_pp_interpol, & ropp_pp_FFT USE star_fsi, ONLY: cos_window IMPLICIT NONE REAL(wp), DIMENSION(:), INTENT(in) :: ang ! central angle between 1st and 2nd radius-vectors REAL(wp), DIMENSION(:,:), INTENT(in) :: full_phase! full phase (m) after propagating to circular orbit[ch, t] REAL(wp), DIMENSION(:,:), INTENT(in) :: amp ! Amplitudes [ch, time] REAL(wp), INTENT(in) :: r1 ! Radius of GNSS (m) REAL(wp), INTENT(in) :: r2 ! Radius of LEO (m) REAL(wp), DIMENSION(:,:), INTENT(inout) :: pp ! Impact parameters (m) REAL(wp), DIMENSION(:,:), INTENT(inout) :: ee ! Bending angles (rad) REAL(wp), DIMENSION(:,:), INTENT(inout) :: aa ! FSI amplitude REAL(wp), INTENT(in) :: pmin ! min impact height REAL(wp), INTENT(in) :: pmax ! max impact height LOGICAL, INTENT(in) :: apply_SH ! flag to apply shadow height determination INTEGER, DIMENSION(:), INTENT(out) :: nout ! array size for impact and bending INTEGER, PARAMETER :: nw = 20 ! COSINE filter window size INTEGER, PARAMETER :: nmax1 = 32768 ! number of data in the interpolated signal INTEGER, PARAMETER :: nmax2 = 65536 ! number of data in the extended signal REAL(wp), PARAMETER :: dsh = 200.0_wp ! Shadow border width (m) CHARACTER(len = 256) :: routine CHARACTER(len = 256) :: file1 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: amp_w ! Amplitudes on new grids [point] REAL(wp), DIMENSION(:), ALLOCATABLE :: ang1 ! open angle grids [point] REAL(wp), DIMENSION(:), ALLOCATABLE :: phs1 ! phase on new grids [point] REAL(wp), DIMENSION(:), ALLOCATABLE :: amp1 ! Amplitudes on new grids [point] REAL(wp), DIMENSION(:), ALLOCATABLE :: phs2 ! pahse on double new grids [point] REAL(wp), DIMENSION(:), ALLOCATABLE :: amp2 ! Amplitudes on double new grids [point] REAL(wp), DIMENSION(:), ALLOCATABLE :: a ! complex signal before FFT [point] REAL(wp), DIMENSION(:), ALLOCATABLE :: u REAL(wp), DIMENSION(:), ALLOCATABLE :: v REAL(wp), DIMENSION(:), ALLOCATABLE :: dop REAL(wp), DIMENSION(:), ALLOCATABLE :: aaa REAL(wp), DIMENSION(:), ALLOCATABLE :: fphs REAL(wp), DIMENSION(:), ALLOCATABLE :: frq REAL(wp), DIMENSION(:), ALLOCATABLE :: tet REAL(wp), DIMENSION(:), ALLOCATABLE :: imp REAL(wp), DIMENSION(:), ALLOCATABLE :: Af REAL(wp), DIMENSION(:), ALLOCATABLE :: PH REAL(wp), DIMENSION(:), ALLOCATABLE :: AH REAL(wp), DIMENSION(:), ALLOCATABLE :: k ! Wave number [channel] REAL(wp) :: w, dx1, dx2, x1, x2, x0 REAL(wp) :: dang, dang1, phsinc, phs21, ang21 REAL(wp) :: teta, teta_start, arg1, arg2, x, y, pha, pp0 INTEGER :: nmax, ic, nc, i, j, kk INTEGER :: i0, i1, ib, nr REAL(wp) :: pminc, Algt, Ashd, Athr, Ascl, dpH ! 1.4 Local arrays CALL message_get_routine(routine) CALL message_set_routine('fsi_leo') ! Determination of data sizes nmax = SIZE(full_phase, 2) ! time dimension nc = SIZE(full_phase, 1) ! 2 frequency channels nr = CEILING(REAL(nmax,wp)/8000.0_wp) ALLOCATE(k(nc)) ! 5.1 Computation of wave vectors k(1) = 2.0_wp*pi*f_L1/C_Light ! wavenubmer in [m-1] k(2) = 2.0_wp*pi*f_L2/C_Light !write(*,*)'nmax, lamda(1), lamda(2), r1, r2',nmax, 2.0_wp*pi/k(1), 2.0_wp*pi/k(2), r1, r2 !write(*,*)'full_phase(:, 1) full_phase(:, nmax)', full_phase(1, 1),full_phase(2, 1), & ! full_phase(1, nmax),full_phase(2, nmax) ! Using COSINE filtering the Amplitude at both edges before FFT ALLOCATE(amp_w(nc, nmax)) dx1 = real(nw, KIND=wp) dx2 = real(nw, KIND=wp) x1 = dx1 x2 = real(2 * nmax, KIND=wp) - dx2 DO i = 1, nmax x0 = real(i, KIND=wp) CALL cos_window(x1, x2, dx1, dx2, x0, w) DO ic=1,nc amp_w(ic, i) = w * amp(ic, i) ENDDO ENDDO ALLOCATE(ang1(nmax1)) ! Discretize the central angle dang=ang(nmax)-ang(1) dang1 = (ang(nmax) - ang(1) ) / (nmax1 - 1) DO i = 1, nmax1 ang1(i) = ang(1) + dang1 * (i - 1) ENDDO Channels: DO ic=1, nc ALLOCATE(phs1(nmax1)) ALLOCATE(amp1(nmax1)) ALLOCATE(phs2(nmax2)) ALLOCATE(amp2(nmax2)) ALLOCATE(a(2 * nmax2)) ALLOCATE(u(-nmax2/2:+nmax2/2)) ALLOCATE(v(-nmax2/2:+nmax2/2)) ALLOCATE(dop(-nmax2/2:+nmax2/2)) ALLOCATE(aaa(-nmax2/2:+nmax2/2)) ALLOCATE(fphs(-nmax2/2:+nmax2/2)) ALLOCATE(frq(-nmax2/2:+nmax2/2)) ALLOCATE(tet(-nmax2/2:+nmax2/2)) ALLOCATE(imp(-nmax2/2:+nmax2/2)) CALL ropp_pp_interpol(ang, ang1, full_phase(ic, :), phs1(:)) CALL ropp_pp_interpol(ang, ang1, amp_w(ic, :), amp1(:)) ! Detrending of the phase phsinc = (phs1(nmax1) - phs1(1) ) / (ang1(nmax1) - ang1(1) ) DO i = 1, nmax1 phs2(i) = phs1(i) - (phs1(1) + dang1 * phsinc * (i - 1) ) amp2(i) = amp1(i) ENDDO !if (ic == 1 ) then ! file1='detrended_phase_L1.dat' ! open(17,file=file1) ! write(17,*)'phsinc =',phsinc ! write(17,*)'dang1 =',dang1 ! write(17,*)' no. angle phase lin_phase detr phase' ! do i=1,nmax1 ! write(17,*)i,ang1(i),phs1(i),dang1*phsinc*(i-1),phs2(i) ! end do ! close(17) !endif !if (ic == 2) then ! file1='detrended_phase_L2.dat' ! open(17,file=file1) ! write(17,*)'phsinc =',phsinc ! write(17,*)'dang1 =',dang1 ! write(17,*)' no. angle phase lin_phase detr phase' ! do i=1,nmax1 ! write(17,*)i,ang1(i),phs1(i),dang1*phsinc*(i-1),phs2(i) ! end do ! close(17) !endif ! Extension of the signal DO i = nmax1 + 1, nmax2 phs2(i) = phs2(nmax1) amp2(i) = 0.0_wp ENDDO ! Calculation of the spectrum of the complex downconverted signal DO i = 1, nmax2 pha = k(ic) * phs2(i) a(2 * i - 1) = amp2(i)*0.1_wp * COS(pha) a(2 * i) = amp2(i)*0.1_wp * SIN(pha) ENDDO CALL ropp_pp_FFT(a, 1) DO i = 1, nmax2 / 2 + 1 j = i - 1 u(j) = a(2 * i - 1) v(j) = a(2 * i) aaa(j) = SQRT(u(j)**2 + v(j)**2) fphs(j) = ATAN2(v(j),u(j)) frq(j) = -1.0_wp*j*pi/dang ENDDO DO i = nmax2 / 2 + 1, nmax2 j = i - nmax2 - 1 u(j) = a(2 * i - 1) v(j) = a(2 * i) aaa(j) = SQRT(u(j)**2 + v(j)**2) fphs(j) = ATAN2(v(j),u(j)) frq(j) = -1.0_wp*j*pi/dang ENDDO ! Calculation of the derivative of the phase of the spectrum (delta_phase = phase(i+1) - phase(i)) DO i = - nmax2 / 2, nmax2 / 2 - 1 x = u(i) * u(i + 1) + v(i) * v(i + 1) y = v(i) * u(i + 1) - u(i) * v(i + 1) dop(i) = ATAN2(-y, x) IF (dop(i) < 0.0_wp) dop(i) = dop(i) + 2.0_wp * pi ENDDO ! 2-point smoothing of doppler DO i = + nmax2 / 2 - 1, - nmax2 / 2 + 1, - 1 dop(i) = (dop(i) + dop(i - 1) ) / 2.0_wp ENDDO dop( - nmax2 / 2) = dop( - nmax2 / 2 + 1) dop( + nmax2 / 2) = dop( + nmax2 / 2 - 1) DO i = + nmax2 / 2 - 1, - nmax2 / 2 + 1, - 1 dop(i) = (dop(i) + dop(i - 1) ) / 2.0_wp ENDDO dop( - nmax2 / 2) = dop( - nmax2 / 2 + 1) dop( + nmax2 / 2) = dop( + nmax2 / 2 - 1) ! Specification of parameters phs21 = nmax2 * (full_phase(ic, nmax) - full_phase(ic, 1) ) / nmax1 ang21 = nmax2 * (ang(nmax) - ang(1) ) / nmax1 !write(*,*)'dang =',dang teta_start=ang(1) DO i= -nmax2/2, +nmax2/2-1 tet(i)=teta_start+dop(i)*ang21/(2.0_wp * pi) imp(i)=(phs21 - i * (2.0_wp * pi)/ k(ic) ) / ang21 ENDDO !IF (ic == 1) THEN ! file1='FFT_phase_amp_dop_freq_L1.txt' !ELSE ! file1='FFT_phase_amp_dop_freq_L2.txt' !ENDIF ! !open(13,file=file1) !write(13,*) r1,r2, pmin, pmax !write(13,*)' N pseu-freq FFT-phase '// & ! ' FFT-dopp theta(rad) impact param FFT ampl' !do i=-nmax2/2, +nmax2/2-1 ! write(13,'(I10,F17.5,3F15.8,2F18.4,2F10.7)') & ! i,frq(i),fphs(i),dop(i),tet(i),imp(i),aaa(i), imp(i)/r1, imp(i)/r2 !end do !close(13) write(*,*)'phs21,ang21',phs21,ang21 !------------------------------------------------------------------------------- ! Determination of shadow border !------------------------------------------------------------------------------- ALLOCATE(Af(nmax2)) ALLOCATE(PH(nmax2)) ALLOCATE(AH(nmax2)) ! make sure the impact heights are monotonous incresing PH(1:nmax2) = imp(+nmax2/2-1:-nmax2/2:-1) AH(1:nmax2) = aaa(+nmax2/2-1:-nmax2/2:-1) i0 = MIN(SUM(MINLOC(PH(:), Pmax > PH(:) .AND. PH(:) > Pmin)), & SUM(MAXLOC(PH(:), Pmax > PH(:) .AND. PH(:) > Pmin))) i1 = MAX(SUM(MINLOC(PH(:), Pmax > PH(:) .AND. PH(:) > Pmin)), & SUM(MAXLOC(PH(:), Pmax > PH(:) .AND. PH(:) > Pmin))) dpH = ABS(PH(i1) - PH(i0))/(i1 - i0) !write( *, *) 'i0, i1, dpH=', i0, i1, dpH IF (apply_SH) THEN ! Determination of light and shadow amplitudes Algt = SQRT(SUM(AH(:)**2, & Mask = (PH(:) > Pmax-5000.0_wp) .AND. (PH(:) < Pmax)) / & COUNT(Mask = (PH(:) > Pmax-5000.0_wp) .AND. (PH(:) < Pmax))) Ashd = SQRT(SUM(AH(:)**2, & Mask = (Pmin-1000.0 < PH(:) .AND. PH(:) < Pmin)) / & COUNT(Mask = (Pmin-1000.0 < PH(:) .AND. PH(:) < Pmin))) ! Determination of threshold and scaling amplitude Athr = 0.5_wp*(Algt + Ashd) Ascl = MIN(Athr, AH(I1) - Ashd) ! Computation of correlation with step function Af(i1) = Ascl DO i=i1-1,i0,-1 Ascl = MIN(Athr, AH(i)) - Ashd Af(i) = Af(i+1) + Ascl ENDDO DO i=i1,i0,-1 Af(i) = Af(i)/SQRT(REAL(i1+1-i)) ENDDO ! Determination of shadow zone border and ! shifting it to nearest point of reduced grid ib = i0 + SUM(MAXLOC(Af(i0:i1))) - 1 ib = ib + NINT(dSh/dpH) ib = MIN(i1-nr, MAX(i0,ib)) ib = nr*CEILING(REAL(ib-1)/REAL(nr)) + 1 ! ib = nr*CEILING(REAL(ib-1)/REAL(nr)) + 5 PminC = PH(ib) !WRITE(*, *) Algt, Ashd, Athr, Ascl, ib, PminC ELSE ! 8.6 Setting min impact height if no automatic shadow zone determination PminC = Pmin ENDIF DEALLOCATE(Af) DEALLOCATE(PH) DEALLOCATE(AH) ! Calculation of the bending angles and impact parameters kk = 1 DO i = - nmax2 / 2, + nmax2 / 2 - 1 arg1 = imp(i) / r1 arg2 = imp(i) / r2 ! IF (imp(i) >= pmin .and. imp(i) <= pmax) THEN IF (imp(i) >= PminC .and. imp(i) <= pmax) THEN IF ( (arg1 > 1.0_wp) .or. (arg2 > 1.0_wp) ) then pp(ic, kk) = 0.0_wp ee(ic, kk) = 0.0_wp aa(ic, kk) = 0.0_wp ELSE pp(ic, kk) = imp(i) ee(ic, kk) = ASIN(arg1) + ASIN(arg2) + tet(i) - pi aa(ic, kk) = aaa(i) ENDIF kk = kk + 1 ENDIF ENDDO nout(ic) = kk - 1 write(*,*) 'kk:', kk !DO i=1,nmax2,40 ! write(*,*) dop(i) * ang21 / (2.0_wp * pi) !ENDDO DEALLOCATE(phs1) DEALLOCATE(amp1) DEALLOCATE(phs2) DEALLOCATE(amp2) DEALLOCATE(a) DEALLOCATE(u) DEALLOCATE(v) DEALLOCATE(dop) DEALLOCATE(aaa) DEALLOCATE(fphs) DEALLOCATE(frq) DEALLOCATE(tet) DEALLOCATE(imp) ENDDO Channels DEALLOCATE(k) DEALLOCATE(ang1) END SUBROUTINE fsi_leo