program f90_RANDWAVE
    implicit none
    integer::iii,i,n
    real(4)::dt,a,b,c,accmax
    real(4),allocatable::y1(:),y2(:),y3(:)
    real(4),allocatable::z1(:),z2(:),z3(:)
    character::fnameW*50,linebuff*300

    dt=0.01
    accmax=100.0

    do iii=1,3
        select case(iii)
            case(1)
                a=5.0;b=15.0;c=30.0 ;fnameW='inp_ran1.csv'
                n=int(c/dt);allocate(y1(1:n));allocate(z1(1:n))
                call WAVE(n,dt,a,b,c,y1,z1)
            case(2)
                a=5.0;b=25.0;c=60.0 ;fnameW='inp_ran2.csv'
                n=int(c/dt);allocate(y2(1:n));allocate(z2(1:n))
                call WAVE(n,dt,a,b,c,y2,z2)
            case(3)
                a=5.0;b=35.0;c=120.0;fnameW='inp_ran3.csv'
                n=int(c/dt);allocate(y3(1:n));allocate(z3(1:n))
                call WAVE(n,dt,a,b,c,y3,z3)
        end select
        open (12, file=fnameW,status='replace')
        write(12,'(a)') trim(fnameW)
        write(linebuff,*) 'dt,',dt
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
        write(linebuff,*) 'ndata,',n+1
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
        write(12,'(a)') '0,0'
        do i=1,n
            if(iii==1)write(linebuff,'(f10.3,",",f10.6)') y1(i)*accmax,z1(i)
            if(iii==2)write(linebuff,'(f10.3,",",f10.6)') y2(i)*accmax,z2(i)
            if(iii==3)write(linebuff,'(f10.3,",",f10.6)') y3(i)*accmax,z3(i)
            call del_spaces(linebuff)
            write(12,'(a)') trim(linebuff)
        end do
        close(12)
    end do

    stop

contains

    subroutine WAVE(n,dt,a,b,c,x,z)
        !wave multiplied by the envelope function for -1.0<y<1.0
        integer,intent(in)::n
        real(4),intent(in)::dt,a,b,c
        real(4),intent(out)::x(1:n),z(1:n)

        real(4)::alpha,coef,tt

        call random_number(x)
        alpha=-log(0.1)/(c-b)
        do i=1,n
            tt=dt*real(i)
            if(tt<a)coef=(tt/a)**2.0
            if(a<=tt.and.tt<b)coef=1.0
            if(b<=tt)coef=exp(-alpha*(tt-b))
            x(i)=(x(i)-0.5)*2.0*coef
            z(i)=coef
        end do
    end subroutine WAVE


    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_RANDWAVE