program rstf c ////////////////////////////////////////////////////////////////////////////////// c c The code for determination of relative source time function by c Non Negative Least Square (NNLS) inversion techniques. c c Code autor: V. Plicka, vladimir.plicka@mff.cuni.cz c Date: May 2019 c Version: v.10.3 c Changes compared to the previous version: c adding smoothing according JZ, time_fixed_smooth.for c c the inversion is NOT done for whole seismograms. c Prescribed time window is applayed to seismograms. c w_cut1 ... left cutting bondary c w_cut2 ... right cutting bondary c cut_length ... duration of cosine smmothing window. c The cosine filter is applayed before w_cut1 and after w_cut2 c c Removing computing VR for individual time intervals. This is not important in case of c dealing with time window cutting c c ////////////////////////////////////////////////////////////////////////////////// parameter (nsmpls = 16384, nsmpls3 = 49152, nsrch = 400) c parameter (nsmpls = 25000, nsmpls3 = 49152, nsrch = 400) complex ai, exp2 complex wx_f(nsmpls), wy_f(nsmpls), wz_f(nsmpls), & ux_f(nsmpls), uy_f(nsmpls), uz_f(nsmpls), & hx(nsmpls), hy(nsmpls), hz(nsmpls), & hx_simple(nsmpls), hy_simple(nsmpls), hz_simple(nsmpls), & rf_mat(nsmpls) real shift, T, dt, pi, pi2, Mo_main, Mo_egf, & tw1, tw2, win_left, win_right real uxr(nsmpls), uyr(nsmpls), uzr(nsmpls), & ux(nsmpls), uy(nsmpls), uz(nsmpls), & x_synt(nsmpls), y_synt(nsmpls), z_synt(nsmpls), & x_synt_simple(nsmpls), y_synt_simple(nsmpls), & z_synt_simple(nsmpls), & wxr(nsmpls), wyr(nsmpls), wzr(nsmpls), & wx(nsmpls), wy(nsmpls), wz(nsmpls), & d(nsmpls), d_simple(nsmpls) real x(nsmpls, nsrch), y(nsmpls, nsrch), z(nsmpls, nsrch), & rst_mat(nsmpls) integer np, nn, i, j, k, cmp(3), & RMo_choise, keydis double precision rnorm, a(nsmpls3,nsrch), b(nsmpls3), & ampl(nsrch), ampl_simple(nsrch) c declaration needed for XAPIIR, filtration subroutine c------------------------------------------------ CHARACTER*2 TYPE, APROTO INTEGER NSAMPS, PASSES, IORD REAL*4 TRBNDW, AS, FLO, FHI, TS, SN(30), SD(30) LOGICAL ZP c------------------------------------------------ open(11, file='main.dat') open(12, file='egf.dat') open(13, file="input.inp") open(14, file='stf.dat') open(15, file='main_synt.dat') open(20, file='main_disp.dat') open(21, file='egf_disp.dat') open(22, file='main_fil.dat') open(23, file='egf_fil.dat') open(25, file='sp_main.dat') open(26, file='sp_egf.dat') open(28, file='sp_main_synt.dat') open(29, file='sp_ratio.dat') open(31, file="vr_final.dat") c open(32, file="vr_complete.dat") open(111,file='inv.dat',status='unknown') open(112,file='inv_cumsum.dat',status='unknown') open(113,file='mo2mo.dat',status='unknown') open(114,file='win_limits.dat',status='unknown') read(13,*) ! Station code, not used by program, for plotting only (rstf.gpl) read(13,*)np_origin, dt ! Number of time samples, time step in seconds read(13,*)w_cut1, w_cut2, cut_length ! Window cut in seconds. The seismograms out othe window are 0. read(13,*)cmp(1),cmp(2),cmp(3) ! 1/0 = Y/No for EW, NS, Z components read(13,*)t_start, t_end ! Time limits (Start, End) for SEARCH read(13,*)win_left, win_right ! In seconds, must be inside . Forward simulation is done for this interval. Outside this window, the weights are=0 read(13,*)shift ! Shift of triangles read(13,*)Mo_main, Mo_egf ! Mo of mainshock and aftershock read(13,*)RMo_choise ! 1/0 Yes/No moment ratio read(13,*)weigh_choise ! weigh_choise, 1/0 Yes/No read(13,*)keydis ! 1/0 Yes/No displacement/velocity, cumsum process, INTEGER read(13,*)keyfil ! keyfil, 0=BP band pass, 1=LP low pass, 2=HP high pass read(13,*)f1, f4 ! Low, high cut in Hertz c read(13,*) c read(13,*)tw1, tw2 ! Definition of the time windows for computing VR, (o-s)^2 and o^2 read(13,*)keysmo ! Smoothing 1/0 Yes/No, Keysmo, ai = (0.,1.) pi = 3.141592 pi2 = 2. * pi T = 1/f4 c------------------------------------------------ RMo = Mo_main/Mo_egf nc = 3 ! Number of components c------------------------------------------------ c Is required, if displacement, computed by cumsum. uxxx=0. uxxy=0. uxxz=0. wxxx=0. wxxy=0. wxxz=0. c------------------------------------------------ c------------------------------------------------ c Set up max amplitudes to zero. Weighting amax_x = 0. amax_y = 0. amax_z = 0. c------------------------------------------------ c Reading the mainshock=U1 and egf=U0. do i = 1, np_origin read(11,*) time, uxr(i), uyr(i), uzr(i) ! main read(12,*) time, wxr(i), wyr(i), wzr(i) ! egf enddo do i = 1, np_origin tt = float(i-1)*dt if (abs(uxr(i)) .ge. amax_x) amax_x = abs(uxr(i)) if (abs(uyr(i)) .ge. amax_y) amax_y = abs(uyr(i)) if (abs(uzr(i)) .ge. amax_z) amax_z = abs(uzr(i)) if(keydis.eq.1) then ! displacement is required, cumsum uxxx = uxxx +uxr(i)*dt uxxy = uxxy +uyr(i)*dt uxxz = uxxz +uzr(i)*dt uxr(i)=uxxx uyr(i)=uxxy uzr(i)=uxxz wxxx = wxxx +wxr(i)*dt wxxy = wxxy +wyr(i)*dt wxxz = wxxz +wzr(i)*dt wxr(i)=wxxx wyr(i)=wxxy wzr(i)=wxxz write(20,*)tt,uxr(i),uyr(i),uzr(i) ! main_dis.dat write(21,*)tt,wxr(i),wyr(i),wzr(i) ! egf_disp.dat endif enddo close(11) ! main.dat close(12) ! egf.dat close(13) ! input.dat close(20) ! main_dis.dat close(21) ! egf_dis.dat write(*,*) c------------------------------------------------ c Weighing YES/NO c Found the maximum from amax_x(y,z) if (weigh_choise.eq.1) then write(*,*) "Weighing" xyz_max = amax_x if (amax_y .ge. xyz_max) xyz_max = amax_y if (amax_z .ge. xyz_max) xyz_max = amax_z weight_x = xyz_max / amax_x weight_y = xyz_max / amax_y weight_z = xyz_max / amax_z else weight_x = 1 weight_y = 1 weight_z = 1 endif c------------------------------------------------ if(keyfil.eq.3) goto 40 c XAPIIR implementation c------------------------------------------------ c NSAMPS=np c TS=dt APROTO='BU' ! Butterworth TRBNDW=0. ! not used for 'BU' AS=0. ! not used for 'BU' IORD=4 ! number of poles (4-5 recommended) the 4 for compatibility with ISOLA GUI Butt amplitude response 4th, power 8th PASSES=1 ! (1 for forward; 2 for forward and reverse = zero phase) MAX_NT=8192 if(keyfil.eq.0) TYPE='BP' ! band pass if(keyfil.eq.1) TYPE='LP' ! low pass if(keyfil.eq.2) TYPE='HP' ! high pass FLO=f1 ! ignored for LP if(FLO.eq.0.) FLO=1./(dt*8192.) ! BandPass cannot start at 0 FHI=f4 call XAPIIR(uxr, np_origin, APROTO, TRBNDW, AS, IORD, + TYPE, FLO, FHI, dt, PASSES, MAX_NT) call XAPIIR(uyr, np_origin, APROTO, TRBNDW, AS, IORD, + TYPE, FLO, FHI, dt, PASSES, MAX_NT) call XAPIIR(uzr, np_origin, APROTO, TRBNDW, AS, IORD, + TYPE, FLO, FHI, dt, PASSES, MAX_NT) call XAPIIR(wxr, np_origin, APROTO, TRBNDW, AS, IORD, + TYPE, FLO, FHI, dt, PASSES, MAX_NT) call XAPIIR(wyr, np_origin, APROTO, TRBNDW, AS, IORD, + TYPE, FLO, FHI, dt, PASSES, MAX_NT) call XAPIIR(wzr, np_origin, APROTO, TRBNDW, AS, IORD, + TYPE, FLO, FHI, dt, PASSES, MAX_NT) c------------------------------------------------ 40 write(*,*)"" c------------------ S T A R T ------------------------------ c Select time window for inversion c Cosine tapperis applayed at the edge of the time window w_cut0 = w_cut1 - cut_length w_cut3 = w_cut2 + cut_length iw_cut0 = w_cut0/dt+1.1 iw_cut1 = w_cut1/dt+1.1 iw_cut2 = w_cut2/dt+1.1 iw_cut3 = w_cut3/dt+1.1 ncut = cut_length/dt dwind = pi/ncut c------------------------------------------------ c Some checking must be included !!!!!!!!!!!! if (w_cut0 .lt. 0.) then write(*,*)"Left cutting parameter is less than 0." write(*,*)"Increase w_cut1, or decrease cut_length" stop endif if (w_cut3 .gt. np_origin) then write(*,*)"Right cutting parameter is greater than np_origin." write(*,*)"Decrease w_cut2, or decrease cut_length" stop endif ki = 0 ! re-indexed the time series arrays according the cut window do i = 1, np_origin c re-indexed the time series arrays according the cut window if ( (i.ge.iw_cut0) .and. (i.le.iw_cut3) ) then ki = ki + 1 ux(ki) = uxr(i) uy(ki) = uyr(i) uz(ki) = uzr(i) wx(ki) = wxr(i) wy(ki) = wyr(i) wz(ki) = wzr(i) endif enddo np = ki c New limits for cosine time smoothing parameters iw_new0 = 1 iw_new1 = iw_cut1 - iw_cut0 iw_new2 = iw_cut2 - iw_cut0 iw_new3 = np c------------------------------------------------ c Applying time window cut on seismograms. The pre defined time window is inverted only, c no whole seismograms. do i = 1, np ux(i) = ux(i) * & smooth(i, iw_new0, iw_new1, iw_new2, iw_new3, dwind) uy(i) = uy(i) * & smooth(i, iw_new0, iw_new1, iw_new2, iw_new3, dwind) uz(i) = uz(i) * & smooth(i, iw_new0, iw_new1, iw_new2, iw_new3, dwind) wx(i) = wx(i) * & smooth(i, iw_new0, iw_new1, iw_new2, iw_new3, dwind) wy(i) = wy(i) * & smooth(i, iw_new0, iw_new1, iw_new2, iw_new3, dwind) wz(i) = wz(i) * & smooth(i, iw_new0, iw_new1, iw_new2, iw_new3, dwind) enddo c------------------ E N D ------------------------------ c write(*,*) np_origin, iw_cut0, iw_cut3 do i = 1, np_origin time = float(i-1) * dt if ( (i.ge.iw_cut0) .and. (i.le.iw_cut3) ) then k = i - iw_cut0 + 1 write(22,*) time, ux(k),uy(k),uz(k) write(23,*) time, wx(k),wy(k),wz(k) else write(22,*) time, 0., 0., 0. write(23,*) time, 0., 0., 0. endif enddo close(22) ! main_fil.dat close(23) ! egf_fil.dat do i = 1, np wx_f(i) = cmplx(wx(i),0.) wy_f(i) = cmplx(wy(i),0.) wz_f(i) = cmplx(wz(i),0.) ux_f(i) = cmplx(ux(i),0.) uy_f(i) = cmplx(uy(i),0.) uz_f(i) = cmplx(uz(i),0.) c write(22,*) (i-1)*dt, ux(i),uy(i),uz(i) c write(23,*) (i-1)*dt, wx(i),wy(i),wz(i) enddo c Move to frequency domain c------------------------------------------------ ni = aint(alog(float(np))/alog(2.))+1 nt = 2**ni nt2 = nt/2 ntm = nt2+1 ntt = nt+2 df = 1./(dt*float(nt)) np1 = np+1 do i = np1, nt wx_f(i) = cmplx(0., 0.) wy_f(i) = cmplx(0., 0.) wz_f(i) = cmplx(0., 0.) ux_f(i) = cmplx(0., 0.) uy_f(i) = cmplx(0., 0.) uz_f(i) = cmplx(0., 0.) enddo call FCOOLR(NI, wx_f, -1.) call FCOOLR(NI, wy_f, -1.) call FCOOLR(NI, wz_f, -1.) call FCOOLR(NI, ux_f, -1.) call FCOOLR(NI, uy_f, -1.) call FCOOLR(NI, uz_f, -1.) c Spectrum of main shock and EGF c------------------------------------------------ do i=1,ntm freq = float(i-1)*df sp_ux = dt*cabs(ux_f(i)) sp_uy = dt*cabs(uy_f(i)) sp_uz = dt*cabs(uz_f(i)) sp_wx = dt*cabs(wx_f(i)) sp_wy = dt*cabs(wy_f(i)) sp_wz = dt*cabs(wz_f(i)) if (wx_f(i) .eq. (0.,0.) ) then rat_x = dt*cabs(ux_f(i)/1e-8) else rat_x = dt*cabs(ux_f(i)/wx_f(i)) endif if (wy_f(i) .eq. (0.,0.) ) then rat_y = dt*cabs(uy_f(i)/1e-8) else rat_y = dt*cabs(uy_f(i)/wy_f(i)) endif if (wz_f(i) .eq. (0.,0.) ) then rat_z = dt*cabs(uz_f(i)/1e-8) else rat_z = dt*cabs(uz_f(i)/wz_f(i)) endif write(25,'(4(1x,e12.6))')freq, sp_ux, sp_uy, sp_uz write(26,'(4(1x,e12.6))')freq, sp_wx, sp_wy, sp_wz write(29,'(4(1x,e12.6))')freq, rat_x, rat_y, rat_z enddo close(25) ! sp_main.dat close(26) ! sp_egf.dat close(29) ! sp_ratio.dat c------------------------------------------------ c The convolution between EGF and shifted triangles c------------------------------------------------ t_wind = t_end - t_start nn = int(t_wind/shift + 1.1) if (nn.gt.nsrch) then write(*,*) "Number of triangles exceed dimension limit" stop endif write(*,*)"Number of searched triagles= ", nn do j = 1, nn ! nn pocet trojuhelniku = pocet neznamych tau = t_start + float(j-1) * shift c write(*,*)tau do i = 1,ntm fr = float(i-1) * df exp1 = pi2 * fr * tau exp2 = -ai * cmplx(exp1, 0.) c hx(i) = wx_f(i) * d_f(i) * exp(exp2) c hy(i) = wy_f(i) * d_f(i) * exp(exp2) c hz(i) = wz_f(i) * d_f(i) * exp(exp2) hx(i) = wx_f(i) * exp(exp2) hy(i) = wy_f(i) * exp(exp2) hz(i) = wz_f(i) * exp(exp2) enddo c Transfer to time domain c------------------------------------------------ do i = 2, nt2 k = ntt - i hx(k) = conjg(hx(i)) hy(k) = conjg(hy(i)) hz(k) = conjg(hz(i)) enddo call FCOOLR(ni, hx, +1.) call FCOOLR(ni, hy, +1.) call FCOOLR(ni, hz, +1.) do i = 1, np time = float(i-1) * dt x(i,j) = real(hx(i))/float(nt) y(i,j) = real(hy(i))/float(nt) z(i,j) = real(hz(i))/float(nt) enddo enddo ! Loop end through shiftted triangles c Creation of matrix A from x(i,j), y(i,j) a z(i,j) c------------------------------------------------ inc = 0 do k =1, nc if (cmp(k) .eq. 1) inc = inc + 1 enddo if ( inc.eq.0 ) then write(*,*)"No selected components" stop endif mm = inc*np do j = 1, nn ! Loop over sfifted triangles (unknows) l = 0 do k = 1, nc ! Loop over selected components do i = 1, np ! Loop over time samples if ( (k .eq. 1) .and. (cmp(k).eq. 1) ) then l=l+1 a(l,j) = x(i,j) * weight_x endif if ( (k .eq. 2) .and. (cmp(k).eq. 1) ) then l=l+1 a(l,j) = y(i,j) * weight_y endif if ( (k .eq. 3) .and. (cmp(k).eq. 1) ) then l=l+1 a(l,j) = z(i,j) * weight_z endif enddo enddo enddo if(keysmo.eq.1) then write(*,*) write(*,*) "Enter the smooth constant: <1,100>" write(*,*) "The larger constant means large smoothing." read(*,*) smooth_inver amamax=maxval(a)*smooth_inver endif c------------------------------------------------ c Creation of vector b from ux(i), uy(i) and uz(i) c------------------------------------------------ l = 0 do k = 1, nc do i = 1, np if ( (k .eq. 1) .and. (cmp(k).eq. 1) ) then l = l+1 b(l) = ux(i) * weight_x endif if ( (k .eq. 2) .and. (cmp(k).eq. 1) ) then l = l+1 b(l) = uy(i) * weight_y endif if ( (k .eq. 3) .and. (cmp(k).eq. 1) ) then l = l+1 b(l) = uz(i) * weight_z endif enddo enddo c------------------------------------------------ c------------------------------------------------ c Free or fixed moment ratio if (RMo_choise .eq. 1) then write(*,*) write(*,*)"Fixed moment ratio" write(*,*) write(*,*)"Enter the weighting of moments ratio constraint =" read(*,*) const write(*,*) c c const = 1 mm = mm +1 do j = 1, nn c a(mm,j) = 1. a(mm,j) = const*1./RMo enddo c b(mm) = RMo b(mm) = const* 1. endif c------------------------------------------------ c------------------------------------------------ c Adding Smoothing, May 2019, START j = mm if(keysmo.eq.1) then c j = mm NTR2 = nn ! NTR2 is number of weights (unknown) c second finite difference do irow=2,NTR2-1 j=j+1 b(j)= 0. enddo endif numdat=j c write(*,*)j, numdat, mm c------------------------------------------------ j = mm if(keysmo.eq.1) then c j = mm NTR2 = nn ! NTR2 is number of weights (unknown) c second finite difference do irow=2,NTR2-1 ! adding NTR2-2 rows j=j+1 do jaux=1,NTR2 a(j,jaux)=0. enddo ! MULTIPLICATION by amamax is ESSENTIAL !!!!!!!!! if not, smoothing may overweight inversion a(j,irow+1)=+1 *amamax! * weig(ir,icom) ! *1.e10 a(j,irow )=-2 *amamax! * weig(ir,icom) ! *1.e10 a(j,irow-1)=+1 *amamax! * weig(ir,icom) ! *1.e10 enddo ! over it = ober rows endif numdat=j c write(*,*)j, numdat, mm c Adding Smoothing, May 2019, END c------------------------------------------------ c Calling NNLS c call tinv(a, b, ampl,mode,rnorm,mm,nn) call tinv(a, b, ampl,mode,rnorm,numdat,nn) c CAUTION !!!! tinv is calling NNLS; c write(*,*) 'mode= ',mode write(*,*)"" if(mode.ne.1) then write(*,*) 'mode not 1, problem in nnls' stop endif c------------------------------------------------ c------------------------------------------------ c Saving the results c write(111,*) c write(111,*) 'Inversion result:' c write(111,*) 'moment coefficients of the triangles ampl(i)' c write(111,'(6(1x,e12.6))') (ampl(n),n=1,nn) C Check if sum of ampl(n) is similar to RMo ratio suma_ampl = 0 ampl_max = 0 do n=1,nn if (ampl(n) .gt. ampl_max )ampl_max = ampl(n) suma_ampl = suma_ampl +ampl(n) enddo write(*,*)"Pre-define the Mo/mo ratio =", RMo write(*,*)"Sum of weights, must be similar to Mo/mo = ", suma_ampl write(*,*)"" write(113,*)RMo, suma_ampl, const ! mo2mo.dat write(114,*)win_left, 0, int(ampl_max + 1) ! win_limits.dat write(114,*)win_right, 0, int(ampl_max + 1) ! win_limits.dat c Write normalize cumulative sum of weights cumsum_ampl = 0. do n=1,nn tau = t_start + float(n-1) * shift cumsum_ampl = cumsum_ampl + ampl(n) write(111,*)tau, ampl(n) ! inv.dat write(112,*)tau, cumsum_ampl/suma_ampl ! inv_cumsum.dat enddo close(111) close(112) close(113) close(114) c------------------------------------------------ c------------------------------------------------ c Forward problem c Convolution of the EGF with founded RSTF (weights = ampl(j)) c------------------------------------------------ c Define the time window c The weights inside this window will be used in forward problem. c The weights outside this window will be substitute by 0. n_win_left = int((win_left - t_start)/shift + 1.1) n_win_right = int((win_right - t_start)/shift + 1.1) c write(*,*) n_win_left, n_win_right do j = 1, nn if ( (j .lt. n_win_left) .or. (j .gt. n_win_right) ) then ampl_simple(j) = 0. else ampl_simple(j) = ampl(j) endif enddo do i = 1,ntm hx(i) = (0.,0.) hy(i) = (0.,0.) hz(i) = (0.,0.) hx_simple(i) = (0.,0.) hy_simple(i) = (0.,0.) hz_simple(i) = (0.,0.) fr = float(i-1) * df do j = 1, nn ! nn pocet trojuhelniku = pocet neznamych tau = t_start + float(j-1) * shift c write(*,*)tau exp1 = pi2 * fr * tau exp2 = -ai * cmplx(exp1, 0.) hx(i) = hx(i) + wx_f(i) * ampl(j) * exp(exp2) hy(i) = hy(i) + wy_f(i) * ampl(j) * exp(exp2) hz(i) = hz(i) + wz_f(i) * ampl(j) * exp(exp2) hx_simple(i) = hx_simple(i) + & wx_f(i) * ampl_simple(j) * exp(exp2) hy_simple(i) = hy_simple(i) + & wy_f(i) * ampl_simple(j) * exp(exp2) hz_simple(i) = hz_simple(i) + & wz_f(i) * ampl_simple(j) * exp(exp2) enddo enddo do i =1, ntm freq = float(i-1) * df sp_hx = dt * cabs(hx(i)) sp_hy = dt * cabs(hy(i)) sp_hz = dt * cabs(hz(i)) write(28,'(4(1x,e12.6))')freq, sp_hx, sp_hy, sp_hz ! sp_main.dat enddo do i = 2, nt2 k = ntt - i hx(k) = conjg(hx(i)) hy(k) = conjg(hy(i)) hz(k) = conjg(hz(i)) hx_simple(k) = conjg(hx_simple(i)) hy_simple(k) = conjg(hy_simple(i)) hz_simple(k) = conjg(hz_simple(i)) enddo call FCOOLR(ni, hx, +1.) call FCOOLR(ni, hy, +1.) call FCOOLR(ni, hz, +1.) call FCOOLR(ni, hx_simple, +1.) call FCOOLR(ni, hy_simple, +1.) call FCOOLR(ni, hz_simple, +1.) c------------------------------------------------ c Computing RSTF as a sum of shifted triangles weighted according the ampl(j) c do i = 1, np time = float(i-1) * dt d(i) = 0. d_simple(i) = 0. do j = 1, nn ashift = float(j-1) * shift d(i) = d(i) + ampl(j) * fsource(time, T, ashift) d_simple(i) = d_simple(i) + & ampl_simple(j) * fsource(time, T, ashift) enddo enddo call integ(d, suma, np, dt, nsmpls) call integ(d_simple, suma_simple, np, dt, nsmpls) c write(*,*)kk,suma c write(*,*) "RMo, suma, RMo/suma ",RMo, suma, RMo/suma if (suma.eq.0)then c scale=RMo/1 scale=Mo_main else c scale=RMo/suma scale=Mo_main/suma endif if (suma_simple.eq.0)then scale_simple=Mo_main else scale_simple=Mo_main/suma_simple endif np_wind = (t_end - t_start)/dt +1 nT2 = T/2/dt +1 do i =nT2, np_wind+nT2 time = float(i) * dt d(i) = d(i)*scale c d_simple(i) = d_simple(i)*scale_simple d_simple(i) = d_simple(i)*scale write(14,*)time+t_start-T/2, d(i), d_simple(i) ! stf.dat enddo c------------------------------------------------ c write synthetic main shock c computed VARRED norm1 = 0. ! observed^2 anorm3 = 0. ! (observed - synthetic)^2 anorm1x = 0. anorm1y = 0. anorm1z = 0. anorm3x = 0. anorm3y = 0. anorm3z = 0. anorm3x_s = 0. anorm3y_s = 0. anorm3z_s = 0. do i = 1, np time = float(i-1) * dt x_synt(i) = real(hx(i))/float(nt) y_synt(i) = real(hy(i))/float(nt) z_synt(i) = real(hz(i))/float(nt) x_synt_simple(i) = real(hx_simple(i))/float(nt) y_synt_simple(i) = real(hy_simple(i))/float(nt) z_synt_simple(i) = real(hz_simple(i))/float(nt) anorm1x = anorm1x + ux(i)*ux(i) anorm1y = anorm1y + uy(i)*uy(i) anorm1z = anorm1z + uz(i)*uz(i) anorm3x = anorm3x + (ux(i) - x_synt(i))*(ux(i) - x_synt(i)) anorm3y = anorm3y + (uy(i) - y_synt(i))*(uy(i) - y_synt(i)) anorm3z = anorm3z + (uz(i) - z_synt(i))*(uz(i) - z_synt(i)) anorm3x_s = anorm3x_s & + (ux(i) - x_synt_simple(i))*(ux(i) - x_synt_simple(i)) anorm3y_s = anorm3y_s & + (uy(i) - y_synt_simple(i))*(uy(i) - y_synt_simple(i)) anorm3z_s = anorm3z_s & + (uz(i) - z_synt_simple(i))*(uz(i) - z_synt_simple(i)) anorm1 = anorm1 + ux(i)*ux(i) + uy(i)*uy(i) + uz(i)*uz(i) anorm3 = anorm3 + (ux(i) - x_synt(i))*(ux(i) - x_synt(i)) & + (uy(i) - y_synt(i))*(uy(i) - y_synt(i)) & + (uz(i) - z_synt(i))*(uz(i) - z_synt(i)) anorm3_s = anorm3_s & + (ux(i) - x_synt_simple(i))*(ux(i) - x_synt_simple(i)) & + (uy(i) - y_synt_simple(i))*(uy(i) - y_synt_simple(i)) & + (uz(i) - z_synt_simple(i))*(uz(i) - z_synt_simple(i)) enddo do i = 1, np_origin time = float(i-1) * dt if ( (i.ge.iw_cut0) .and. (i.le.iw_cut3) ) then k = i - iw_cut0 + 1 write(15,'(7(1x,e12.6))')time, & x_synt(k), y_synt(k), z_synt(k), & x_synt_simple(k), y_synt_simple(k), z_synt_simple(k) else write(15,'(7(1x,e12.6))')time, 0., 0., 0., 0., 0., 0. endif enddo close(15) ! main_synt.dat close(28) ! sp_main_synt.dat varredx = 1. - anorm3x/anorm1x varredy = 1. - anorm3y/anorm1y varredz = 1. - anorm3z/anorm1z varred = 1. - anorm3/anorm1 varredx_s = 1. - anorm3x_s/anorm1x varredy_s = 1. - anorm3y_s/anorm1y varredz_s = 1. - anorm3z_s/anorm1z varred_s = 1. - anorm3_s/anorm1 c Computing final VR for used components and individual comp. c------------------------------------------------ write(*,*)"VR for 3-components:", varred write(*,*)"VR for NS:", varredx write(*,*)"VR for EW:", varredy write(*,*)"VR for Z:", varredz write(31,'(4(1x,f7.3))')varredx, varredy, varredz, varred write(31,'(4(1x,f7.3))')varredx_s, varredy_s, varredz_s, varred_s close(31) 812 format('Time wind. ', "VR ", "o^2 ", "(o-s)^2" ) 811 format('o^2 = ', e9.3, ' (o-s)^2 = ', e9.3) c------------------------------------------------ stop end c /////////////////////////////////////////////////////////// c Subroutines & functions c /////////////////////////////////////////////////////////// function fsource(time, tn, shift) real fsource tpul = tn / 2. t0 = shift t1 = shift + tpul t2 = shift + tn a = 2. / tn c a = 1. if ( time .lt. t0 ) fsource = 0. if ( ( time .ge. t0 ) .and. ( time .le. t1 ) ) * fsource = a * ( time - t0 ) / ( t1 - t0 ) if ( ( time .gt. t1 ) .and. ( time .le. t2 ) ) * fsource = a - a * ( time - t1 ) / ( t2 - t1 ) if ( time .gt. (shift+tn) ) fsource = 0. return end subroutine integ(f,suma,np,dt,nsmpls) real f(nsmpls) suma = 0. do i =1, np-1 suma = suma + abs(f(i+1)+f(i))/2. enddo suma = suma * dt return end c------------------------------------------------ c SUBROUTINE FCOOLR(K,D,SN) c FAST FOURIER TRANSFORM OF N = 2**K COMPLEX DATA POINTS c REPARTS HELD IN D(1,3,...2N-1), IMPARTS IN D(2,4,...2N). SUBROUTINE FCOOLR(K,D,SN) DIMENSION INU(20),D(16384) LX=2**K Q1=LX IL=LX SH=SN*6.28318530718/Q1 DO 10 I=1,K IL=IL/2 10 INU(I)=IL NKK=1 DO 40 LA=1,K NCK=NKK NKK=NKK+NKK LCK=LX/NCK L2K=LCK+LCK NW=0 DO 40 ICK=1,NCK FNW=NW AA=SH*FNW W1=COS(AA) W2=SIN(AA) LS=L2K*(ICK-1) DO 20 I=2,LCK,2 J1=I+LS J=J1-1 JH=J+LCK JH1=JH+1 Q1=D(JH)*W1-D(JH1)*W2 Q2=D(JH)*W2+D(JH1)*W1 D(JH)=D(J)-Q1 D(JH1)=D(J1)-Q2 D(J)=D(J)+Q1 20 D(J1)=D(J1)+Q2 DO 29 I=2,K ID=INU(I) IL=ID+ID IF(NW-ID-IL*(NW/IL)) 40,30,30 30 NW=NW-ID 29 CONTINUE 40 NW=NW+ID NW=0 DO 6 J=1,LX IF(NW-J) 8,7,7 7 JJ=NW+NW+1 J1=JJ+1 JH1=J+J JH=JH1-1 Q1=D(JJ) D(JJ)=D(JH) D(JH)=Q1 Q1=D(J1) D(J1)=D(JH1) D(JH1)=Q1 8 DO 9 I=1,K ID=INU(I) IL=ID+ID IF(NW-ID-IL*(NW/IL)) 6,5,5 5 NW=NW-ID 9 CONTINUE 6 NW=NW+ID RETURN END c------------------------------------------------ real function smooth(i, iw_cut0, iw_cut1, iw_cut2, iw_cut3, dwind) integer iw_cut0, iw_cut1, iw_cut2, iw_cut3 if ((i.lt.iw_cut0).or. (i.gt.iw_cut3)) smooth = 0. if ((i.ge.iw_cut1).and.(i.le.iw_cut2)) smooth = 1. if ((i.ge.iw_cut0).and.(i.lt.iw_cut1)) * smooth = 0.5 * (1. - cos(dwind*float(i-iw_cut0))) if ((i.gt.iw_cut2).and.(i.le.iw_cut3)) * smooth = 0.5 * (1. + cos(dwind*float(i-iw_cut2))) return end c------------------------------------------------ include "timefilters.inc"