
	subroutine psv_init

	include 'common.h'

	real np


c source

	a=aa

	if(a.eq.0.)then
	so(1,it0)=1.

	else

	do i=1,nt
	t=(i-it0)*dt

c	so(1,i)=-2.*(a**2)*t*exp(-(a*t)**2)
c	so(1,i)=exp(-(a*t)**2)
c	so(1,i)=(2*(pi*aa*t)**2-1)*exp(-(pi*aa*t)**2)

	np=1.

	if(t.lt.aa.and.t.ge.0)then
            so(1,i) = sin (np*pi*t/aa) - (np/(np+2)) * sin ((np+2)*pi*t/aa)
	endif

	enddo

c normalize
	
	smax=0.
	do i=1,nt
	if(abs(so(1,i)).gt.smax)then
	smax=abs(so(1,i))
	endif
	enddo
	do i=1,nt
	so(1,i)=so(1,i)/smax
	enddo

	if(so(1,1).gt..0001)then
	write(*,*)' Source not small enough at it=1 '
	stop
	endif

	endif


c 

	open(unit=10,file=seisfile(1:lnblnk(seisfile))//'_src')
	do i=1,nt
	write(10,*)so(1,i)
	enddo
	close(10)

c initialize Radius

c	forall(iz=1:nz)vr2(iz)=Rmax-(iz-izfree)*dr
c	forall(iz=1:nz)vr1(iz)=vr2(iz)-.5*dr
c
c	forall(ix=1:nx,iz=1:nz)r1(ix,iz)=vr1(iz)
c	forall(ix=1:nx,iz=1:nz)r2(ix,iz)=vr2(iz)

c initialize theta as a function of source location
c the source location is defined one grid poi nt away from the 
c rotational axis

	forall(ix=1:nx-isx+1)vtheta1(isx-1+ix)=0+(ix-1)*dphi
	forall(ix=1:nx-isx+1)vtheta2(isx-1+ix)=0+(ix-1)*dphi-0.5*dphi
c	forall(ix=1:nx)vtheta1(ix)=0+(ix-1)*dphi
c	forall(ix=1:nx)vtheta2(ix)=0+(ix-1)*dphi-0.5*dphi

c set it anti symmetric around isx-1

	forall(ix=1:isx-1)vtheta1(isx-ix)=-vtheta1(isx+ix)	
	forall(ix=1:isx-1)vtheta2(isx-ix)=-vtheta2(isx+1+ix)	

	forall(ix=1:nx,iz=1:nz)theta1(ix,iz)=vtheta1(ix)
	forall(ix=1:nx,iz=1:nz)theta2(ix,iz)=vtheta2(ix)


c initialization

	u=0.
	r1inv=1./r1
	r2inv=1./r2
	rhoinv=1./rho
	
	where(theta1<>0)
	cot1=cos(theta1)/sin(theta1)
	else where 
	cot1=1.
	end where

	cot2=cos(theta2)/sin(theta2)

c source gauss

	if(alpha<>0.)then

	write(*,*)' alpha = ',alpha
        alpha=1./(alpha*alpha)

	zz0=rmax-sdepth

c        do ix=1,nx
c             do iz=1,nz
c             xx=r2(ix,iz)*sin(theta2(ix,iz))
c             zz=r2(ix,iz)*cos(theta2(ix,iz))-zz0
c		write(*,*)xx,zz
c                gauss(ix,iz)=exp( -alpha *( xx**2 + zz**2))
c             enddo
c        enddo

        forall( ix=1:nx,iz=1:nz)
                gauss(ix,iz)=exp( -alpha *( (r2(ix,iz)*sin(theta2(ix,iz)))**2 
     &                                          +  (r2(ix,iz)*cos(theta2(ix,iz))-zz0)**2))
        end forall



	write(*,*)minval(gauss),' < gauss < ',maxval(gauss)

c	open(unit=10,file=seisfile(1:lnblnk(seisfile))//'_gs')
c	do i=1,nx
c	do j=1,nz
c	write(10,*)gauss(i,j)
c	enddo
c	enddo
c	close(10)

	endif

	
c  absorbing boundary

c     absorb=1.

c         if(nabsorb>0) then
c          forall(ix=1:nabsorb)absorb(ix,:)=ix
c           absorb=exp(-aabsorb*absorb**2 )
c	endif

c Initialize the source distance in x
 


	return
	end
