forward.for

该程序用于两个直立六面体上半无源空间和下半无源空间组合磁异常规则平面网正演.

!     该程序用于两个直立六面体上半无源空间和下半无源空间组合磁异常规则平面网正演.
!     单位:坐标单位m,磁化强度A/m,磁异常nT.
!     观测面坐标:xx,yy,zz;直立六面体坐标x,y,z
      program forward
        real,allocatable::x1(:),y1(:),z1(:),x2(:),y2(:),z2(:),
     &	   xx(:),yy(:),T1(:,:),T2(:,:),T(:,:)
	  integer i,j,nx,ny
	  real Ax,Ay,In,M,xx1,xx2,xxd,yy1,yyd,zz
        character*80 cmdfilename,filename_data,filename_xxyyzz,
     &	  filename_T	   
	  cmdfilename='cmd.txt'

        call read_filename(cmdfilename,filename_data,filename_xxyyzz,
     &	   filename_T)
	  call read_xxyyzz(nx,ny,xx1,xx2,yy1,yy2,zz,filename_xxyyzz)

	  allocate(x1(1:2),y1(1:2),z1(1:2),x2(1:2),y2(1:2),z2(1:2),
     &	       xx(1:nx),yy(1:ny),T1(1:nx,1:ny),T2(1:nx,1:ny),	       
     &           T(1:nx,1:ny))
	  
     	  call read_data(Ax,Ay,In,M,x1,y1,z1,x2,y2,z2,filename_data)  
	  call grid_xxyy(nx,ny,xx1,xx2,yy1,yy2,xx,yy)
	  call calcul_T(Ax,Ay,In,M,x1,y1,z1,xx,yy,zz,nx,ny,T1)
        call calcul_T(Ax,Ay,In,M,x2,y2,z2,xx,yy,zz,nx,ny,T2)
        call combine_T(nx,ny,T1,T2,T)	 
        call output(nx,ny,xx1,xx2,yy1,yy2,xx,yy,T,filename_T)

        deallocate(x1,y1,z1,x2,y2,z2,xx,yy,T1,T2,T)
	end program	
	
!     读入文件名.
	subroutine read_filename(filename,filename_data,filename_xxyyzz,
     &           filename_T)
	  character*(*) filename,filename_data,filename_xxyyzz,filename_T          
	  open(11,file=filename,status='old')
	    read(11,*)filename_data
		read(11,*)filename_xxyyzz
	    read(11,*)filename_T
	  close(11)
	end subroutine    
      
!     读入异常体(直立六面体)的物性参数和空间位置.
      subroutine read_data(Ax,Ay,In,M,x1,y1,z1,x2,y2,z2,filename)
        real e,p,Ax,Ay,In,M,x1(1:2),y1(1:2),z1(1:2),x2(1:2),
     &	   y2(1:2),z2(1:2)
        character*(*)filename
	  e=1e-07
	  open(12,file=filename,status='old')
          read(12,*)Ax,Ay
		read(12,*)In
	    read(12,*)M
	    read(12,*)x1(1),x1(2)
          read(12,*)y1(1),y1(2)
          read(12,*)z1(1),z1(2)
	    read(12,*)x2(1),x2(2)
          read(12,*)y2(1),y2(2)
          read(12,*)z2(1),z2(2)
	    if((Ax+Ay-90)>e.or.(Ax+Ay-270)>e.or.(abs(Ax-Ay)
     &	  -90)>e)write(*,*)'Ax or Ay error'
	  close(12)
	end subroutine

!     读入观测平面坐标范围.
      subroutine read_xxyyzz(nx,ny,xx1,xx2,yy1,yy2,zz,filename)
	  integer nx,ny
        real xx1,xx2,yy1,yy2,zz
	  character*(*)filename
	  open(13,file=filename,status='old')
	    read(13,*)
	    read(13,*)ny,nx
          read(13,*)yy1,yy2
          read(13,*)xx1,xx2
	    read(13,*)zz,zz
        close(13)
	end subroutine

!     将观测面进行规则网格化.
      subroutine  grid_xxyy(nx,ny,xx1,xx2,yy1,yy2,xx,yy)
        integer nx,ny
        real xx1,xxd,yy1,yyd,xx(1:nx),yy(1:ny)
	  xxd=(xx2-xx1)/(nx-1)
	  yyd=(yy2-yy1)/(ny-1)
	  do i=1,nx
	    xx(i)=xx1+(i-1)*xxd
	  end do
	  do i=1,ny
	    yy(i)=yy1+(i-1)*yyd
	  end do
      end subroutine

!     计算异常体在观测面上的磁异常值.
      subroutine calcul_T(Ax,Ay,In,M,x,y,z,xx,yy,zz,nx,ny,T)
	  integer i,j,nx,ny,k,r,s
	  real Ax,Ay,In,M,Mx,My,Mz,zz,Tmin,tx,ty,tz,loge,arct,
     &       x(1:2),y(1:2),z(1:2),xx(1:nx),yy(1:ny),Hax(1:nx,1:ny),
     &	   Hay(1:nx,1:ny),Za(1:nx,1:ny),Vxx(1:nx,1:ny),Vyy(1:nx,1:ny),
     &      Vxy(1:nx,1:ny),Vxz(1:nx,1:ny),Vzz(1:nx,1:ny),Vyz(1:nx,1:ny),
     &      T(1:nx,1:ny)
        tx=cosd(In)*cosd(Ax)
	  ty=cosd(In)*cosd(Ay)  
	  tz=sind(In)
	  Mx=M*tx
	  My=M*ty
	  Mz=M*tz
	  do i=1,nx
	   do j=1,ny
	    Vxx(i,j)=0.0
	    Vyy(i,j)=0.0
	    Vzz(i,j)=0.0
	    Vxy(i,j)=0.0
	    Vxz(i,j)=0.0
	    Vyz(i,j)=0.0
	    do k=1,2
	     do r=1,2
	      do s=1,2
	       Vxx(i,j)=Vxx(i,j)+arct(x(k),y(r),z(s),xx(i),yy(j),zz)
     &		       *(-1.0)**mod(k+r+s,2)
             Vyy(i,j)=Vyy(i,j)+arct(y(r),x(k),z(s),yy(j),xx(i),zz)
     &		       *(-1.0)**mod(k+r+s,2)
             Vzz(i,j)=Vzz(i,j)+arct(z(s),y(r),x(k),zz,yy(j),xx(i))
     &		       *(-1.0)**mod(k+r+s,2)    
             Vxy(i,j)=Vxy(i,j)+loge(x(k),y(r),z(s),xx(i),yy(j),zz)
     &		       *(-1.0)**mod(k+r+s,2)
	       Vxz(i,j)=Vxz(i,j)+loge(x(k),z(s),y(r),xx(i),zz,yy(j))
     &		       *(-1.0)**mod(k+r+s,2)
	       Vyz(i,j)=Vyz(i,j)+loge(z(s),y(r),x(k),zz,yy(j),xx(i))
     &		       *(-1.0)**mod(k+r+s,2)
	       end do
	      end do
	     end do
     	     Hax(i,j)=(Mx*Vxx(i,j)+My*Vxy(i,j)+Mz*Vxz(i,j))
	     Hay(i,j)=(Mx*Vxy(i,j)+My*Vyy(i,j)+Mz*Vyz(i,j))
	     Za(i,j)=(Mx*Vxz(i,j)+My*Vyz(i,j)+Mz*Vzz(i,j))
           T(i,j)=Hax(i,j)*tx+Hay(i,j)*ty+Za(i,j)*tz
	     T(i,j)=T(i,j)*1e+02
	   end do
	  end do
	end subroutine          

!     磁异常自然对数函数.
      function loge(a,b,c,x0,y0,z0)
	  real loge,a,b,c,x0,y0,z0,r0
	  r0=sqrt((a-x0)**2+(b-y0)**2+(c-z0)**2)
	  if(r0+(c-z0)>1e-07)then
	    loge=log(r0+(c-z0))
	  else
	    loge=-log(r0+abs(c-z0))
	  end if
	end function

!     磁异常反正切函数.
      function arct(a,b,c,x0,y0,z0)
	  real pi,arct,a,b,c,x0,y0,z0,r0
	  pi=3.1415926
	  r0=sqrt((a-x0)**2+(b-y0)**2+(c-z0)**2)
	  if(abs(r0*(a-x0))<1e-7)then
	    if(abs((b-y0)*(c-z0))>1e-7)then
	      arct=sign(pi/2,(b-y0)*(c-z0))
	    else
	      arct=0.0
	    end if
	  else
          arct=atan((b-y0)*(c-z0)/(r0*(a-x0)))
	  end if
	  arct=-arct	  
	end function

!     两模型磁异常在观测面的叠加.
      subroutine combine_T(nx,ny,T1,T2,T)
        integer nx,ny
	  real T1(1:nx,1:ny),T2(1:nx,1:ny),T(1:nx,1:ny)
	  do i=1,nx
	    do j=1,ny
	      T(i,j)=T1(i,j)+T2(i,j)
	    end do
	  end do
      end subroutine

!     输出地表各点的组合磁异常值.
      subroutine output(nx,ny,xx1,xx2,yy1,yy2,xx,yy,T,filename)         
	  integer i,j,nx,ny
	  real xx1,xx2,yy1,yy2,Tmin,Tmax,xx(1:nx),yy(1:ny),T(1:nx,1:ny)	  
        character*(*)filename
	  open(14,file=filename)
	  close(14,status='delete') 
	  open(15,file=filename,status='new')
	    do i=1,nx
	      do j=1,ny
	        write(15,*)xx(i),yy(j),T(i,j)
	      end do
	    end do
	  close(22)
 	  Tmin=huge(Tmin)
	  Tmax=-huge(Tmax)
	  do i=1,nx
	    do j=1,ny
	      if(Tmin>T(i,j))Tmin=T(i,j)
            if(Tmax<T(i,j))Tmax=T(i,j)
	    end do
	  end do
	  open(16,file='data.grd')
	  close(16,status='delete') 
	  open(17,file='data.grd',status='new')        
	    write(17,"(a)")"DSAA"
	    write(17,*)ny,nx
	    write(17,*)yy1,yy2
	    write(17,*)xx1,xx2
	    write(17,*)Tmin,Tmax
          do i=1,nx
    	      write(17,*)(T(i,j),j=1,ny)
	    end do
        close(17)          
 	end subroutine