c This programs is used to calculate the sensitivity kernels at a period.
c the input file is sensinp,   run as " sensitivity < sensinp " 

c For each period,the input parameters in sensinp are 
c   (1) average phase velocity of your study area, 
c   (2) spectral file for the singe-frency seismograms windowed and cutted with 
c       the same window length as that applied to real data at this period 
     
c   (3) outfn which stores the sensitivity kernels 
c   (4) smooth length which is used to smooth the sensivitiy kernels and should be 
c       the same length as that used to get phase velocity map from inverted phase velocity at 
c       each grid.  



	parameter( maxnfreq = 30,radius = 6371,
     1             maxix = 400,maxiy = 400 , niter = 11 )

	real    phsens(maxix,maxiy),      ampsens(maxix,maxiy)
	real avgphsens(maxix,maxiy),   avgampsens(maxix,maxiy)
	real wgttemp(maxix,maxiy)
        real scaazi(361)

	real freq(maxnfreq),amplitude(maxnfreq)
	real kk
	character *70 outfn1,outfn2,spcfn,scaazifn


	write(*,*) 'input velocity average phase velocity'
	read(*,*) phvel
	write(*,*) 'please input the spectral file'
	read(*,'(a)') spcfn
	write(*,*) 'type outputfile of the smoothed sensitiviy kernels'
	read(*,*) outfn2
	write(*,*) 'please input smooth len (km) '
	read(*,*)   scalelen

         write(*,*) phvel, spcfn,outfn2

	open(10,file = spcfn)
	open(30,file = outfn2)





	xbeg = -1500
	xend =  1500
	ybeg = -1500
	yend =  1500
	dx   =  10.
	dy   =  10.
	
	nx = (xend-xbeg)/dx + 1
	ny = (yend-ybeg)/dy + 1	 
	
	pi = 4.* atan(1.)   

	open(10,file = spcfn)
	  sumamp = 0.

	  ifreq = 1
 20	  read(10,*,end=40) freq(ifreq),amplitude(ifreq)
           write(*,*)  freq(ifreq),amplitude(ifreq)
	  sumamp = sumamp + amplitude(ifreq)
	   ifreq = ifreq + 1
	   goto 20
 40	   continue
	   nfreq = ifreq - 1
          write(*,*) nfreq

	do ifreq = 1,nfreq
	   amplitude(ifreq) = amplitude(ifreq)/sumamp
	enddo	
	
	do ix = 1,maxix
	do iy = 1,maxiy
	   phsens(ix,iy) = 0.
	avgphsens(ix,iy) = 0.
	   ampsens(ix,iy) = 0. 
	avgampsens(ix,iy) = 0.
	enddo
	enddo

	do ifreq = 1, nfreq
	   write(*,*) ifreq
	 period = 1./freq(ifreq)
	 lamda = phvel*period
         kk = 2*pi/lamda*radius
  	 do ix = 1,nx
	     x = (ix-1)*dx + xbeg
  	     delta1 = x
	   do iy = 1,ny
	     y = (iy-1)*dy + ybeg
	     if(x.eq.0 .and. y.eq.0) then 
	       delta2 = sqrt(dx**2+dy**2)
	     else
	       delta2 = sqrt(x**2+y**2)
	     endif

	    angle = atan2(-y,-x)*180./pi
	    if ( angle .lt. 0) angle = angle + 360.

           iangle = int(angle)


        phsens(ix,iy) = phsens(ix,iy)   + 
     1  amplitude(ifreq)*(-2)*kk**2*sin(kk*(delta1+delta2)/radius+pi/4)
     1                 /sqrt(8*pi*kk*abs(sin(delta2/radius)))
     1                 * ((dx*dy)/radius**2)

	ampsens(ix,iy) = ampsens(ix,iy) +
     1  amplitude(ifreq)*(-2)*kk**2*cos(kk*(delta1+delta2)/radius+pi/4)
     1                 /sqrt(8*pi*kk*abs(sin(delta2/radius)))
     1              * ((dx*dy)/radius**2)
  
                    
	   enddo
	 enddo
	enddo	     




        alpha = 1./( scalelen**2 )

	nxlimit = 20. 	
	nylimit = 20. 
	write(30,*) 301, -1500., 10.
        write(30,*) 301, -1500., 10. 

         do ix = 1, nx
	 do iy = 1, ny
	   x = (ix-1)*dx + xbeg
	   y = (iy-1)*dy + ybeg
            wgtsum = 0
   
	   ixxbeg = ix-nxlimit
	   ixxend = ix+nxlimit
	   iyybeg = iy-nylimit
	   iyyend = iy+nylimit

           if( ix .le. nxlimit    ) ixxbeg = 1
	   if( ix .ge. (nx-nxlimit)) ixxend = nx
           if( iy .le. nylimit    ) iyybeg = 1
	   if( iy .ge. (ny-nylimit)) iyyend = ny
	
           do ixx = ixxbeg,ixxend
	   do iyy = iyybeg,iyyend
	      xx = (ixx-1)*dx + xbeg
	      yy = (iyy-1)*dy + ybeg
              distsq = alpha*((xx-x)**2+(yy-y)**2)
             if( distsq.lt. 80. ) then
               wgttemp(ixx,iyy) = exp(-distsq)
               wgtsum = wgtsum + wgttemp(ixx,iyy)
             else
               wgttemp(ixx,iyy) = 0.0
             endif
           enddo
	   enddo


           do ixx = ixxbeg,ixxend
	      do iyy = iyybeg,iyyend
	       avgphsens(ix,iy)   = avgphsens(ix,iy) + 
     1                     phsens(ixx,iyy)*wgttemp(ixx,iyy)/wgtsum
	      avgampsens(ix,iy)  = avgampsens(ix,iy) + 
     1                    ampsens(ixx,iyy)*wgttemp(ixx,iyy)/wgtsum
              enddo
	   enddo        

	   write(30,*) x,y,avgphsens(ix,iy),avgampsens(ix,iy)
	 enddo
	 enddo		
	 
	  close(10)
	  close(30)


	end