module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
end module defpi


program f90_wave_make
    !*****************************************************************
    ! artificial earthquake motion
    !*****************************************************************
    use defpi
    implicit none
    integer::ndata,i,k,nn
    real(8)::accmax,velmax,dismax
    real(8),allocatable::accinp(:)       !Acceleration fo original wave
    real(8),allocatable::acc(:)          !Acceleration
    real(8),allocatable::vel(:)          !Velocity
    real(8),allocatable::dis(:)          !Displacement
    real(8),allocatable::resaccinp(:)    !Response of acceleration of original wave
    real(8),allocatable::resacc(:)       !Response of acceleration
    real(8),allocatable::resvel(:)       !Response of velocity
    real(8),allocatable::resdis(:)       !Response of displacement
    real(8),allocatable::xr(:)           !Complex Fourier series real part
    real(8),allocatable::xi(:)           !Complex Fourier series imaginary part
    real(8),allocatable::frq(:),Tk(:)    !frequency,period
    real(8),allocatable::rk(:)           !spectrum ratio (target wave/calculated wave)
    real(8)::err1,err2,eps
    real(8)::sum,tt,x1,y1,x2,y2
    integer::iii,maxita                  !Limit value of iteration
    real(8)::hh                          !Damping factor
    real(8)::dt                          !Time increment
    integer::nts                         !Number of data of target spectrum
    real(8),allocatable::perts(:)        !Period of target spectrum
    real(8),allocatable::accts(:)        !Acceleration of target spectrum (function of period)
    real(8),allocatable::tsacc(:)        !Acceleration of target spectrum (function of frequency)
    character::fnameR1*50,fnameR2*50,fnameW1*50,fnameW2*50
    character::dummy*50,linebuff*300

    maxita=1000 !Set limit value of iteration
    hh=0.05     !Set damping factor
    eps=0.05    !Set tolerance

    call getarg(1,fnameR1)  !Input of original wave (acceleration)
    call getarg(2,fnameR2)  !Input of target spectrum
    call getarg(3,fnameW1)  !output of calculated wave
    call getarg(4,fnameW2)  !Output of response spectrum

    open(11,file=fnameR1,status='old')
    read(11,'(a)') linebuff
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    nn=2  !set number of data to powers of 2
    do
        nn=nn*2
        if(ndata*1<=nn)exit
    end do
    allocate(accinp(1:nn))
    allocate(acc(1:nn))
    allocate(vel(1:nn))
    allocate(dis(1:nn))
    allocate(xr(1:nn))
    allocate(xi(1:nn))
    allocate(rk(1:nn/2+1))
    allocate(frq(1:nn/2+1))
    allocate(Tk(1:nn/2+1))
    allocate(resacc(1:nn/2+1))
    allocate(resvel(1:nn/2+1))
    allocate(resdis(1:nn/2+1))
    allocate(tsacc(1:nn/2+1))
    allocate(resaccinp(1:nn/2+1))

    do k=1,nn/2+1 !set frequency
        frq(k)=1.0D0/dt/dble(nn)*dble(k-1)
    end do
    frq(1)=0.5D0*frq(2)
    do k=1,nn/2+1
        Tk(k)=1.0D0/frq(k)
    end do
    do i=1,nn
        accinp(i)=0.0D0
    end do
    do i=1,ndata
        read(11,*) accinp(i)
    end do
    close(11)

    !Input of target spectrun
    open(11,file=fnameR2,status='old')
    read(11,'(a)') linebuff
    read(11,*) dummy,nts
    read(11,*) dummy,dummy
    allocate(perts(1:nts))
    allocate(accts(1:nts))
    do i=1,nts
        read(11,*) perts(i),accts(i)
    end do
    close(11)
    !target acceleration at calculated frequency
    do k=1,nn/2+1 !frequency
        tt=1.0D0/frq(k)
        if(tt<perts(1))then
            x1=log10(perts(1));y1=log10(accts(1))
            x2=log10(perts(2));y2=log10(accts(2))
        end if
        if(perts(nts-1)<tt)then
            x1=log10(perts(nts-1));y1=log10(accts(nts-1))
            x2=log10(perts(nts));y2=log10(accts(nts))
        end if
        if(perts(1)<=tt.and.tt<=perts(nts-2))then
            do i=1,nts-1
                if(perts(i)<=tt.and.tt<=perts(i+1))then
                    x1=log10(perts(i))  ;y1=log10(accts(i))
                    x2=log10(perts(i+1));y2=log10(accts(i+1))
                end if
            end do
        end if
        tsacc(k)=10.0D0**CALY(log10(tt),x1,y1,x2,y2)
    end do

    do i=1,nn
        acc(i)=accinp(i)
    end do
    err1=1.0D0
    !***********************
    do iii=1,maxita
    !***********************
        do i=1,nn
            vel(i)=0.0D0
            dis(i)=0.0D0
        end do
        call IACC(nn,dt,acc,vel,dis,accmax,velmax,dismax)         !Integral of acceleration
        call CRAC(nn,dt,acc,vel,dis,accmax,velmax,dismax)         !Amendment of base line
        call CALRES(nn,nn/2+1,dt,hh,acc,Tk,resacc,resvel,resdis) !calculation of response
        if(iii==1)then
            do k=1,nn/2+1
                resaccinp(k)=resacc(k)
            end do
        end if
        do k=1,nn/2+1
            rk(k)=tsacc(k)/resacc(k)
        end do
        sum=0.0D0
        do k=1,nn/2+1
            sum=sum+(1.0D0-rk(k))**2.0D0
        end do
        err2=sqrt(sum/dble(nn/2+1))
        !***********************
        if(err2<=eps)exit
        if(10<iii.and.err1<err2)exit
        print "('iii=',i0,' error=',e15.7)",iii,err2
        err1=err2
        !***********************
        do i=1,nn
            xr(i)=acc(i)
            xi(i)=0.0D0
        end do
        call FFT(nn,xr,xi)  !Fourier transform
        do i=1,nn  !To divide return value by number of data
            xr(i)=xr(i)/dble(nn)
            xi(i)=xi(i)/dble(nn)
        end do
        do k=1,nn/2+1 !Update of Fourier coefficient
            xr(k)=xr(k)*rk(k)
            xi(k)=xi(k)*rk(k)
        end do
        do k=nn/2+2,nn !To set Fourier coefficient to all degree
            xr(k)=xr(nn-k+2)
            xi(k)=-xi(nn-k+2)
        end do
        xi(1)=0.0D0       !To ensure the characteristic of imaginary part
        xi(nn/2+1)=0.0D0  !To ensure the characteristic of imaginary part
        xr(1)=0.0D0       !To set average value of wave to zero
        do i=1,nn  !To inverse the sign of imaginary part
            xr(i)=xr(i)
            xi(i)=-xi(i)
        end do
        Call FFT(nn,xr,xi)  !Inverse Fourier transform
        do i=1,nn
            acc(i)=0.0D0
        end do
        x1=0.0
        do i=1,ndata
            acc(i)=xr(i)
            if(x1<abs(acc(i)))x1=abs(acc(i)) !To search max.acceleration
        end do
        do i=1,ndata !Amendment of max.acceleration
            acc(i)=acc(i)/x1*tsacc(nn/2+1)
        end do
        acc(1)=0.0D0   !Zero set of initial value
    !***********************
    end do
    !***********************


    open (12, file=fnameW1, status='replace')
    linebuff='acceleration'
    write(12,'(a)') trim(linebuff)
    write(linebuff,*) 'dt,',dt
    call del_spaces(linebuff)
    write(12,'(a)') trim(linebuff)
    write(linebuff,*) 'nn,',nn
    call del_spaces(linebuff)
    write(12,'(a)') trim(linebuff)
    do i=1,nn
        write(linebuff,*) accinp(i),',',acc(i)
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
    end do
    close(12)

    open (12, file=fnameW2, status='replace')
    write(dummy,'(f6.3)') err2
    linebuff='spectrum(error='//trim(dummy)//')'
    call del_spaces(linebuff)
    write(12,'(a)') trim(linebuff)
    write(linebuff,*) 'nn/2=,',nn/2
    call del_spaces(linebuff)
    write(12,'(a)') trim(linebuff)
    write(12,'(a)') 'frq(Hz),period(sec),org(gal),obtained(gal),target(gal)'
    do k=1,nn/2+1
        write(linebuff,*) frq(k),',',1.0D0/frq(k),',',resaccinp(k),',',resacc(k),',',tsacc(k)
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
    end do
    close(12)

    stop

contains

    real(8) function CALY(xx,x1,y1,x2,y2)
        !**************************************************
        !Simple regression : yy=a*xx+b
        !**************************************************
        real(8),intent(in)::xx,x1,y1,x2,y2
        real(8)::a,b

        a=(y2-y1)/(x2-x1)
        b=(x2*y1-x1*y2)/(x2-x1)
        CALY=a*xx+b
    end function CALY


    subroutine FFT(nn,xr,xi)
    !**************************************************
    !Fast Fourier Transform, Inverse transform
    !**************************************************
    !nn   : Number of data (powers of 2)
    !xr() : Real part of i/o data
    !xi() : Imaginary part of i/o data
        integer,intent(in)::nn
        real(8),intent(inout)::xr(1:nn)
        real(8),intent(inout)::xi(1:nn)

        integer::g,h,i,j,k
        integer::l,m,n,p,q
        real(8)::a,b,xd
        real(8),allocatable::s(:)
        real(8),allocatable::c(:)

        n=nn
        allocate(s(1:n/2+1))
        allocate(c(1:n/2+1))
        i=0;j=0;k=0;l=0;p=0;h=0;g=0;q=0
        m=int(log(dble(n))/log(2.0D0)+1.0D0)
        a=0.0D0
        b=pi*2.0D0/dble(n)
        do i=0,n/2
            s(i+1)=sin(a)
            c(i+1)=cos(a)
            a=a+b
        end do
        l=n
        h=1
        do g=1,m
            l=l/2
            k=0
            do q=1,h
                p=0
                do i=k,l+k-1
                    j=i+l
                    a=xr(i+1)-xr(j+1)
                    b=xi(i+1)-xi(j+1)
                    xr(i+1)=xr(i+1)+xr(j+1)
                    xi(i+1)=xi(i+1)+xi(j+1)
                    if(p==0)then
                        xr(j+1)=a
                        xi(j+1)=b
                    else
                        xr(j+1)=a*c(p+1)+b*s(p+1)
                        xi(j+1)=b*c(p+1)-a*s(p+1)
                    end if
                    p=p+h
                end do
                k=k+l+l
            end do
            h=h+h
        end do
        j=n/2
        do i=1,n-1
            k=n
            if(j<i)then
                xd=xr(i+1)
                xr(i+1)=xr(j+1)
                xr(j+1)=xd
                xd=xi(i+1)
                xi(i+1)=xi(j+1)
                xi(j+1)=xd
            end if
            k=k/2
            do while(j>=k)
                j=j-k
                k=k/2
                if(k==0)exit
            end do
            j=j+k
        end do
    end subroutine FFT


    subroutine NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,dt,ff)
        !*****************************************************************
        !Newmark Beta Method
        !*****************************************************************
        real(8),intent(out)::acc2,vel2,dis2
        real(8),intent(in)::acc1,vel1,dis1
        real(8),intent(in)::am,ac,ak,dt,ff
!        real(8)::beta=1.0D0/4.0D0  !Average acceleration method
        real(8)::beta=1.0D0/6.0D0  !Linear acceleration method

        acc2=(ff-ac*(vel1+0.5D0*dt*acc1)-ak*(dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1))&
             /(am+ac*0.5D0*dt+ak*beta*dt*dt)
        vel2=vel1+0.5D0*dt*(acc1+acc2)
        dis2=dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1+beta*dt*dt*acc2
    end subroutine NEWMARK


    subroutine CALRES(ndata,nnn,dt,hh,acc,Tk,resacc,resvel,resdis)
        !*****************************************************************
        !Calculation of response
        !*****************************************************************
        integer,intent(in)::ndata,nnn
        real(8),intent(in)::dt,hh
        real(8),intent(in)::acc(1:ndata),Tk(1:nnn)
        real(8),intent(out)::resacc(1:nnn),resvel(1:nnn),resdis(1:nnn)

        integer::i,j,k,md
        real(8)::acc2,vel2,dis2
        real(8)::acc1,vel1,dis1
        real(8)::am,ac,ak,f1,f2,fm,ddt
        real(8)::racc,rvel,rdis
        real(8)::acmax,vemax,dimax

        md=5
        ddt=dt/dble(md)
        do k=1,nnn
            am=1.0D0                       !Mass
            ak=4.0D0*pi*pi*am/Tk(k)/Tk(k)  !Stiffness
            ac=2.0D0*hh*sqrt(ak*am)        !Damping constant

            acc1=0.0D0
            vel1=0.0D0
            dis1=0.0D0
            acmax=acc1
            vemax=vel1
            dimax=dis1
            do i=1,ndata
                f1=0.0D0
                if(1<=i)f1=-am*acc(i-1)
                f2=-am*acc(i)
                do j=1,md
                    fm=f1+(f2-f1)/dble(md)*dble(j)
                    call NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,ddt,fm)
                    acc1=acc2
                    vel1=vel2
                    dis1=dis2
                end do
                racc=acc2+acc(i)     !Absolute acceleration
                rvel=vel2            !Relative velocity
                rdis=dis2            !Relative displacement
                if(acmax<abs(racc))acmax=abs(racc)
                if(vemax<abs(rvel))vemax=abs(rvel)
                if(dimax<abs(rdis))dimax=abs(rdis)
            end do
            resacc(k)=acmax
            resvel(k)=vemax
            resdis(k)=dimax
        end do
    end subroutine CALRES


    subroutine IACC(nn,dt,acc,vel,dis,accmax,velmax,dismax)
        !*****************************************************************
        !Integral of acceleration
        !*****************************************************************
        !nn    : Number of time history data of acceleration
        !dt    : Time increment                    (input data)
        !acc() : Time history data of Acceleration (input data)
        !vel() : Time history data of Velocity     (output data)
        !dis() : Time history data of Displacement (output data)
        !accmax: Maximum value of acceleration     (output data)
        !velmax: Maximum value of velocity         (output data)
        !dismax: Maximum value of displacement     (output data)

        integer,intent(in)::nn
        real(8),intent(in)::dt
        real(8),intent(in)::acc(1:nn)
        real(8),intent(out)::vel(1:nn),dis(1:nn)
        real(8),intent(out)::accmax,velmax,dismax

        integer::i

        accmax=0.0
        velmax=0.0
        dismax=0.0
        vel(1)=0.0
        dis(1)=0.0
        do i=2,nn
            vel(i)=vel(i-1)+(acc(i-1)+acc(i))*dt/2.0D0
            dis(i)=dis(i-1)+vel(i-1)*dt+(acc(i-1)/3.0D0+acc(i)/6.0D0)*dt*dt
            if(accmax<abs(acc(i)))accmax=abs(acc(i))
            if(velmax<abs(vel(i)))velmax=abs(vel(i))
            if(dismax<abs(dis(i)))dismax=abs(dis(i))
        end do
    end subroutine IACC


    subroutine CRAC(nn,dt,acc,vel,dis,accmax,velmax,dismax)
        !*****************************************************************
        !Amendment of base line
        !*****************************************************************
        !nn    : Number of time history data of acceleration
        !dt    : Time increment                    (input data)
        !acc() : Time histori data of acceleration (input & output data)
        !vel() : Time history of velocity          (input data)
        !dis() : Time history data of displacement (input & output data)
        !accmax: Maximum value of acceleration     (input data)
        !velmax: Maximum value of velocity         (input data)
        !dismax: Maximum value of displacement     (input data)

        integer,intent(in)::nn
        real(8),intent(in)::dt,accmax,velmax,dismax
        real(8),intent(in)::vel(1:nn)
        real(8),intent(inout)::dis(1:nn)
        real(8),intent(inout)::acc(1:nn)

        integer::i
        real(8)::tt,t,sum,a0,a1,acmax,coef

        tt=dble(nn-1)*dt
        t=0.0D0
        do i=1,nn
            dis(i)=dis(i)*(3.0D0*tt-2.0D0*t)*t*t
            t=t+dt
        end do
        sum=(dis(1)+dis(nn))/2.0D0
        do i=2,nn-1
            sum=sum+dis(i)
        end do
        sum=sum*dt
        a1=28.0D0/13.0D0/tt/tt*(2.0D0*vel(nn)-15.0D0/tt**5.0D0*sum)
        a0=vel(nn)/tt-a1/2.0D0*tt
        t=0.0D0
        acmax=0.0D0
        do i=1,nn
            acc(i)=acc(i)-a0-a1*t
            if(acmax<abs(acc(i)))acmax=abs(acc(i))
            t=t+dt
        end do
        coef=accmax/acmax
        do i=1,nn
            acc(i)=acc(i)*coef
        end do
    end subroutine CRAC


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


end program f90_wave_make

