c  program timedomain
c     time-domain radiated pulse structure
      parameter(nt=16384,nt1=8192,nnz=5)
      real s(nt),s1(nt),rw,c0,freq1,freqn,dz
      complex gr(nt1)
      integer nos,N,np,nzw,nf
      character*20 dtg
      call datime(dtg)
      write(*,'(1X,''START: '',A/)') dtg
      open(unit=1,status='old',file='pulse.in')
      open(unit=2,status='old',file='swampgreen.dat',access='direct',
     &RECL=8,form='binary')
      open(unit=3,status='replace',file='pulse.out',
     &form='binary',RECL=4)
      open(unit=4,status='replace',file='test.out')
c      open(unit=5,status='old',file='ramfreq.in')
      call pulse(nt,s,freq,np,nos,N,nr,nz,nrw,rw,c0,nf,
     &dz,freq1,freqn)
      isign=(1)
      call four1(s,N,isign)
c      call filter(nt,s,freq,np,nos,N,nl,nh)
      do 40 i=1,nz,1
      do 41 jj=1,2*N,1
      s1(jj)=s(jj)
  41  continue
      nzw=i
      call propag(freq,np,nt,nt1,s1,gr,N,nr,nz,nrw,nzw,
     &nf,freq1,freqn,rw,c0)
      isign=(-1)
      call four1(s1,N,isign)
      if(nzw.eq.nnz) then
      jj=0
       do 65 j=1,2*N,2
      amp=sqrt(s1(j)*s1(j)+s1(j+1)*s1(j+1))
      write(4,*)rw/c0+1./freq/np*jj,amp/N
      jj=jj+1
  65  continue
      end if
      call out(nt,s1,N,nzw,freq,np,rw,c0,dz)
 40   continue
      call datime(dtg)
      write(*,'(1X,''STOP: '',A/)') dtg
      close(1)
      close(2)
      close(3)
      close(4)
      close(5)
      end

      subroutine pulse(nt,s,freq,np,nos,N,nr,nz,nrw,rw,c0,
     &nf,dz,freq1,freqn)
      real s(nt),freq,dt,t,rw,c0,dz,freqr,freq1,freqn
      integer nf, iii
cc    Modified by UTSAV/KILARI  
      print*, 'subroutine pulse'
cc      read(1,*)
      read(1,*)freq,np,nos,N,nr,nz,nrw,rw,c0,dz
      freq1=freq-freq/nos
      freqn=freq+freq/nos
      nf=0
      do 3333 iii=0,N/2,1
      freqr=iii*freq*np/N
      if(freqr.lt.freq1) goto 3333
      if(freqr.gt.freqn) goto 3333
      nf=nf+1
 3333 continue     
      ii=0
      pi=4.0*atan(1.0)
      dt=1/freq/np
      do 10 i=1,2*N,2
      t=dt*ii
      ts=1/freq*nos
      if(t.le.ts) then
      s(i)=cos(2*pi*freq*t)
      s(i+1)=0.0
      else
      s(i)=0.0
      s(i+1)=0.0
      end if
      ii=ii+1
 10   continue
      return
      end
      
      subroutine filter(nt,s,freq,np,nos,N,nl,nh)
      real s(nt),buf
      integer nos,N,np
      nl=ifix((1.0-1.0/nos)*N/np)
      nh=ifix((1.0+1.0/nos)*N/np)
      do 20 i=1,N+1,2
      if(i.lt.nl*2) then
      s(i)=0.0
      s(i+1)=0.0
      s(2*N-i+1)=0.0
      s(2*N-i)=0.0
      else if (i.gt.nh*2) then
      s(i)=0.0
      s(i+1)=0.0
      s(2*N-i+1)=0.0
      s(2*N-i)=0.0
      end if
 20   continue
      return
      end

      subroutine propag(freq,np,nt,nt1,s1,gr,N,nr,nz,nrw,nzw,
     &nf,freq1,freqn,rw,c0)
      real s1(nt),freq,freqn,freq1,dfreq,
     &t0,freqc,rw,c0
      complex gr(nt1),cbuf,cbuf1,cbuf2
      integer*4 nrec
      integer j,np,nn,nf
      nn=nf
      t0=rw/c0
      do 29 j=1,nt1,1
      gr(j)=cmplx(0.0,0.0)
  29  continue
      do 30 j=1,nn,1
      nrec=nr*nz*(j-1)+(nrw-1)*nz+nzw
      read(2,rec=nrec) gr(j)
 30   continue
c      if(nzw.eq.80) then
c      do 65 j=1,nn,1
c      amp=sqrt(s(j)*s(j)+s(j+1)*s(j+1))
c       amp=s(j)
c      amp=cabs(gr(j))
c      write(4,*)345.0+460.*3./2048*(j-1),amp
c 65   continue
c     end if
      ii=1
      ii1=0
      pi=4.0*atan(1.0)
      dt=1/freq/np
      do 50 i=1,N+1,2
      freqc=1.0/N/dt*ii1
      if(freqc.lt.freq1) then
      s1(i)=0.0
      s1(i+1)=0.0
      s1(2*N-i+1)=0.0
      s1(2*N-i)=0.0
      else if (freqc.gt.freqn) then
      s1(i)=0.0
      s1(i+1)=0.0
      s1(2*N-i+1)=0.0
      s1(2*N-i)=0.0
      else
      cbuf=cmplx(s1(i),s1(i+1))
      cbuf1=cbuf*gr(ii)
c      cbuf=cmplx(0.,t0*2*pi*freqc)
c      cbuf2=cexp(cbuf)
c      cbuf1=cbuf1*cbuf2
      s1(i)=real(cbuf1)
      s1(2*N-i)=real(cbuf1)
c       s1(2*N-i)=0.0
      s1(i+1)=imag(cbuf1)
      s1(2*N-i+1)=(-1.0)*imag(cbuf1)
c      s1(2*N-i+1)=0.0
      ii=ii+1
      end if
      ii1=ii1+1
 50   continue
      return
      end

      subroutine out(nt,s1,N,nzw,freq,np,rw,c0,dz)
      integer nzw,np
      real s1(nt),freq,rw,c0,dz
      jj=0
      do 60 j=1,2*N,2
      amp=(sqrt(s1(j)*s1(j)+s1(j+1)*s1(j+1)))/N
      amp1=(20.)*log10(amp+1.e-20)
      write(3)(-1.0)*(nzw-1)*dz,rw/c0+1./freq/np*jj,amp1
c      write(3)(-1.0)*(nzw-1)*dz,rw/c0+1./freq/np*jj,amp
      jj=jj+1
  60  continue
      return
      end 

cc      subroutine datime(dtg)
cc      character date*11, months(12)*3, time*11
cc      integer*2 day,hour,minute,second,i100th,month,year
cc      data months / 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
cc     &'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec' /
cc      call getdat (year,month,day)
cc      call gettim (hour,minute,second,i100th)
cc      write(date,'(I2,1X,A3,1X,I4)') day,months(month),year
cc      write(time,'(2(I2,1H:),I2)') hour,minute,second
cc      write(dtg,'(A11,1X,A8)') date,time
cc      return
cc      end

cc old routine datime replaced by this one which
cc is obtained after mofifying datime in swamp.f 
      subroutine datime(dtg)
      character dtg*(*)
      character date*11, months(12)*3, time*11
      integer*2 day,hour,minute,second,i100th,month,year
      data months / 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
     &'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec' /
cc    ADDED date_time() instead of getdat() BY UTSAV/KILARI FOR SOLARIS
	integer date_time(8)
	character*10 b(3)
 	character newdate*11
	character newtime*8
	integer newmonth
cc       call date (newdate) 
cc      call gettim (hour,minute,second,i100th)
cc      write(date,'(I2,1X,A3,1X,I4)') day,months(month),year
	call date_and_time(b(1),b(2),b(3),date_time)
	newmonth = date_time(2)
      write(newdate,'(I2,1X,A3,1X,I4)') date_time(3),months(newmonth),date_time(1) 
      write(newtime,'(2(I2,1H:),I2)') date_time(5),date_time(6),date_time(7)
cc      write(time,'(2(I2,1H:),I2)') hour,minute,second
cc      write(dtg,'(A11,1X,A8)') date,time
      write(dtg,'(A11,1X,A8)') newdate,newtime
      return
      end