continuation.for

该程序用于规则平面磁异常延拓

!     该程序用于规则平面网磁异常延拓.
      program continuation
        real,allocatable::x(:),y(:),w(:,:),T(:,:),Timag(:,:)
	  parameter(pi=3.1415926)
	  integer i,j,n,m,n0,n1,n2,n3,m0,m1,m2,m3
	  real z,z1,xmin,xmax,ymin,ymax,Tmin,Tmax
        character*80 cmdfilename,filename_data,filename_in,filename_out  
	  cmdfilename='cmd.txt'

        call read_filename(cmdfilename,filename_data,filename_in,
     &	   filename_out)
        call read_in(n,m,xmin,xmax,ymin,ymax,z,z1,filename_in)
        call calculate(n,n0,n1,n2,n3)
	  call calculate(m,m0,m1,m2,m3)

        allocate(x(n0:n3),y(m0:m3),w(n0:n3,m0:m3),T(n0:n3,m0:m3),
     &	       Timag(n0:n3,m0:m3))
        Timag=0.0

	  call read_data(n0,n1,n2,n3,m0,m1,m2,m3,x,y,T,filename_data)
        call exp_cos(n0,n1,n2,n3,m0,m1,m2,m3,T)  
        call fft2(T,Timag,n3-n0+1,m3-m0+1,2)        
        call frequency(pi,n0,n1,n2,n3,m0,m1,m2,m3,x,y,w)
	  call continue_h(n0,n3,m0,m3,z,z1,w,T,Timag)
	  call fft2(T,Timag,n3-n0+1,m3-m0+1,1)
	  call minmax_T(n0,n1,n2,n3,m0,m1,m2,m3,T,Tmin,Tmax)
	  call output(n,m,n0,n1,n2,n3,m0,m1,m2,m3,x,y,T,
     &       xmin,xmax,ymin,ymax,Tmin,Tmax,filename_out)
     	  
        deallocate(x,y,w,T,Timag)
	end program	
	
!     读入文件名.
	subroutine read_filename(filename,filename_data,filename_in,
     &	       filename_out)
	  character*(*)filename,filename_data,filename_in,filename_out	       
	  open(11,file=filename,status='old')
	    read(11,*)filename_data
		read(11,*)filename_in
	    read(11,*)filename_out
	  close(11)
	end subroutine
	
!     读入观测数据个数n,m.
      subroutine read_in(n,m,xmin,xmax,ymin,ymax,z,z1,filename)
	  character*(*)filename
	  integer n,m
 	  real xmin,xmax,ymin,ymax,z
	  open(12,file=filename,status='old')
	    read(12,*)
	    read(12,*)m,n
	    read(12,*)ymin,ymax
	    read(12,*)xmin,xmax
	    read(12,*)z,z
	    read(12,*)z1
	end subroutine

!	计算扩边数据个数n0,n1,n2,n3(m0,m1,m2,m3).
      subroutine calculate(n,n0,n1,n2,n3)
	  integer n,n0,n1,n2,n3,k
        k=int(log10(1.0*n)/log10(2.0)+0.5)+1
	  n0=0
	  n1=(2**k-n)/2
	  n2=n1+n-1
	  n3=2**k-1
	end subroutine
	
!     输入磁异常观测数据.
      subroutine read_data(n0,n1,n2,n3,m0,m1,m2,m3,x,y,T,filename)
	  integer n0,n1,n2,n3,m0,m1,m2,m3
	  real x(n0:n3),y(m0:m3),T(n0:n3,m0:m3)
	  character*(*) filename
	  open(13,file=filename,status='old')
	    do i=n1,n2
	      do j=m1,m2
	        read(13,*)x(i),y(j),T(i,j)
	       end do
	    end do
	  close(13)	
	end subroutine
	    
!     对观测数据进行扩边处理.
      subroutine exp_cos(n0,n1,n2,n3,m0,m1,m2,m3,T)
	  integer n0,n1,n2,n3,m0,m1,m2,m3
	  real delt,T(n0:n3,m0:m3),Tmin,Tmax
        delt=(T(n1,m1)+T(n2,m1)+T(n1,m2)+T(n2,m2))/4.0
	  do j=m0,m3
	    T(n0,j)=delt
	    T(n3,j)=delt
	  end do
	  do i=n0,n3
	    T(i,m0)=delt
	    T(i,m3)=delt
	  end do
        do j=m1,m2
	    do i=n0+1,n1-1
	      T(i,j)=T(n0,j)+cosd(90.0*(i-n1)/(n0-n1))*(T(n1,j)-T(n0,j))
	    end do
	    do i=n2+1,n3-1
	      T(i,j)=T(n3,j)+cosd(90.0*(i-n2)/(n3-n2))*(T(n2,j)-T(n3,j))
          end do
        end do
	  do i=n0+1,n3-1
	    do j=m0+1,m1-1
	      T(i,j)=T(i,m0)+cosd(90.0*(j-m1)/(m0-m1))*(T(i,m1)-T(i,m0))
	    end do
	    do j=m2+1,m3-1
	      T(i,j)=T(i,m3)+cosd(90.0*(j-m2)/(m3-m2))*(T(i,m2)-T(i,m3))
	    end do
        end do
	end subroutine

!     二维FFT变换.
      subroutine fft2(datareal,dataimag,m,n,nf)
        real datareal(m,n),dataimag(m,n)
        real,allocatable::treal(:),timag(:)
        nmmax=max(m,n)

        allocate(treal(nmmax),timag(nmmax),stat=ierr)
	  if(ierr/=0)write(*,*)'数组分配错误,stop!'

        do i=1,m
          if(n/=1)then
            do j=1,n
              treal(j)=datareal(i,j)
              timag(j)=dataimag(i,j)
            end do
            call fft(treal,timag,n,nf)
            do j=1,n
              datareal(i,j)=treal(j)
              dataimag(i,j)=timag(j)
            end do
          end if
	  end do
        do j=1,n
          if(m/=1)then
            do i=1,m
              treal(i)=datareal(i,j)
              timag(i)=dataimag(i,j)
            end do
            call fft(treal,timag,m,nf)
            do i=1,m
              datareal(i,j)=treal(i)
              dataimag(i,j)=timag(i)
            end do
          end if
	  end do

        deallocate(treal,timag,stat=ierr)
        
      end subroutine

!     一维FFT变换.
	subroutine fft(xreal,ximag,n,nf)
        real xreal(n),ximag(n)
        nu=int(log(float(n))/0.693147+0.001)
	  n2=n/2
	  nu1=nu-1
	  f=float((-1)**nf)
	  k=0
        DO l=1,nu
          DO while (k<n)
            do i=1,n2
              p=ibitr(k/2**nu1,nu)
              arg=6.2831853*p*f/float(n)
              c=cos(arg)
              s=sin(arg)
              k1=k+1
              k1n2=k1+n2
              treal=xreal(k1n2)*c+ximag(k1n2)*s
              timag=ximag(k1n2)*c-xreal(k1n2)*s
              xreal(k1n2)=xreal(k1)-treal
              ximag(k1n2)=ximag(k1)-timag
              xreal(k1)=xreal(k1)+treal
              ximag(k1)=ximag(k1)+timag
              k=k+1
            end do
            k=k+n2
          end do
          k=0
          nu1=nu1-1
          n2=n2/2
        end do
        do k=1,n
          i=ibitr(k-1,nu)+1
          if(i>k)then
            treal=xreal(k)
            timag=ximag(k)
            xreal(k)=xreal(i)
            ximag(k)=ximag(i)
            xreal(i)=treal
            ximag(i)=timag
          end if
        end do
        if(nf/=1)then
          do i=1,n
            xreal(i)=xreal(i)/float(n)
            ximag(i)=ximag(i)/float(n)
          end do
        end if
      end subroutine

!     外部函数.
      function ibitr(j,nu)
	  j1=j
	  itt=0
        do i=1,nu
          j2=j1/2
          itt=itt*2+(j1-2*j2)
          ibitr=itt
          j1=j2
        end do
      end function

!     计算频率U,V,W.
      subroutine frequency(pi,n0,n1,n2,n3,m0,m1,m2,m3,x,y,w)
	  integer i,j,n,m,n0,n1,n2,n3,m0,m1,m2,m3
	  real pi,delx,dely,x(n0:n3),y(m0:m3),
     &	   u(n0:n3),v(m0:m3),w(n0:n3,m0:m3)
	  n=n3-n0+1
	  m=m3-m0+1
	  delx=(x(n2)-x(n1))/(n2-n1)/2
	  dely=(y(m2)-y(m1))/(m2-m1)/2
        do i=0,n/2
	    u(i)=i/(n*delx)*2*pi
	  end do
	  do i=n/2+1,n-1
	    u(i)=(i-n)/(n*delx)*2*pi
        end do
	  do j=0,m/2
	    v(j)=j/(m*dely)*2*pi
	  end do
	  do j=m/2+1,m-1
	    v(j)=(j-m)/(m*dely)*2*pi
	  end do
	  do i=0,n-1
	    do j=0,m-1
	      w(i,j)=sqrt(u(i)**2+v(j)**2)*2*pi
	    end do
	  end do
      end subroutine

!     计算位场延拓(延拓因子G与频谱相乘).
      subroutine continue_h(n0,n3,m0,m3,z,z1,w,T,Timag)
	  integer n0,n3,m0,m3
	  real(8) G
	  real z,z1,Tmin,Tmax,w(n0:n3,m0:m3),T(n0:n3,m0:m3),
     &	   Timag(n0:n3,m0:m3)
	  do i=n0,n3
	    do j=m0,m3
	      G=exp(-1.0*w(i,j)*(z-z1))
	      T(i,j)=T(i,j)*G
		  Timag(i,j)=Timag(i,j)*G	      
	    end do
        end do
	end subroutine

!     计算磁异常的最大值和最小值.
      subroutine minmax_T(n0,n1,n2,n3,m0,m1,m2,m3,T,Tmin,Tmax)
	  integer n0,n1,n2,n3,m0,m1,m2,m3
	  real Tmin,Tmax,T(n0:n3,m0:m3)
	  Tmin=huge(Tmin)
	  Tmax=-huge(Tmax)
	  do i=n1,n2
	    do j=m1,m2
	      if(Tmin>T(i,j))Tmin=T(i,j)
            if(Tmax<T(i,j))Tmax=T(i,j)
	    end do
	  end do
	end subroutine

!     输出延拓后的磁异常值.
      subroutine output(n,m,n0,n1,n2,n3,m0,m1,m2,m3,x,y,T,
     &           xmin,xmax,ymin,ymax,Tmin,Tmax,filename)
	  integer n,m,n0,n1,n2,n3,m0,m1,m2,m3
	  real Tmin,Tmax,x(n0:n3),y(m0:m3),T(n0:n3,m0:m3)
        character*(*)filename
	  open(14,file=filename)
	  close(14,status='delete') 
	  open(15,file=filename,status='new') 
	    write(15,"(a)")"DSAA"
	    write(15,*)m,n
	    write(15,*)ymin,ymax
	    write(15,*)xmin,xmax    
	    write(15,*)Tmin,Tmax
          do i=n1,n2
    	      write(15,*)(T(i,j),j=m1,m2)
	    end do
        close(15)
 	end subroutine