! 瞬变电磁一维正演
        subroutine Tem_calculate(p,n,t,u,n1,a,i0)
        implicit none
        integer  i,j,n,n1,mm
        real(8) drt,pi,w,a,i0,u0,ff
        real(8) p(2*n-1),u(n1),wsl(-149:100),wl(140)
        real t(n1)
        complex(8) hw
        pi=3.1415926
        u0=4.e-7*pi
        mm=140
        open(120,file='sin250.txt')
        do i=-149,100
        read(120,*)wsl(i)
        enddo
        close(120)
        print*,wsl
        open(130,file='J1_coef.txt')   !J1系数
        do i=1,mm
        read(130,*)wl(i)
        enddo
        close(130)
        print*,wl
        drt=alog(10.0)/20.0            !alog(x)是单精度自然对数 dlog(x)是双精度自然对数 clog(x)是复数度自然对数 log(x)表示自然对数
        u=0
        ff=u0/4.0/pi*(2*pi*i0*a*a*u0/5.0)**(2.0/3.0)!场值
        do i=1,n1
           do j=-149,100,1
              w=exp(j*drt)/t(i)
              call frequency(hw,w,wl,mm,p,n,i0,a)
              u(i)=u(i)+imag(hw)*wsl(j)
           enddo   
           u(i)=u(i)/t(i)*sqrt(pi/2)
           u(i)=abs(u(i)*2/pi*u0)
           u(i)=ff/t(i)*(1.0/t(i)/u(i))**(2.0/3.0)!衰减电压值
        enddo
        print*,u
        pause
       endsubroutine
       !***********************************************
       !频率域一维正演
        subroutine frequency(hw,w,wl,mm,p,n,I0,a)
        implicit none
        external ths
        integer n,mm,i,j    
        real(8) u0,pi,s,I0,a,a0,w
        complex(8)  z(n),zs(n),hw,k(n),u1(n),aa,ths,z0
        real(8)  r(mm),wl(mm),p(2*n-1)
        pi=3.1415926
        u0=4.e-7*pi
        a0=-7.91001919
        s=8.7967143957e-2
        do i=1,mm
           r(i)=10**(a0+s*(i-1))/a
        enddo
         hw=0.0
       do j=1,mm
         k(n)=cmplx(0,-1*w*u0/p(n))!p(i)即第i层电阻率 p(i+n)即第i层厚度 complex(a,b)为转换为a+bi的复数
         z0=cmplx(0,-1*w*u0/r(j))
         u1(n)=sqrt(r(j)**2+k(n))
         z(n)=cmplx(0,-1*w*u0)/u1(n)
         zs(n)=z(n)
         do i=n-1,1,-1
           k(i)=cmplx(0,-1*w/p(i)*u0)
           u1(i)=sqrt(r(j)**2+k(i))
           z(i)=cmplx(0,-1*w*u0)/u1(i)
           aa=p(i+n)*u1(i)
           zs(i)=z(i)*(zs(i+1)+z(i)*ths(aa))/(z(i)+zs(i+1)*ths(aa)) 
         enddo
          hw=hw+r(j)*(zs(1)/(zs(1)+z0)-0.5)*wl(j)
       enddo
       hw=hw*i0
        endsubroutine
        !***********************************************************
        function ths(aa)
        IMPLICIT NONE
        complex(8) ths,aa
        ths=exp(-2*aa)
        ths=(1-ths)/(1+ths)
        endfunction