论文要求：

1、         研究一维地电模型正演程序代码；（1.得出结果需要哪些数据和文件？）

2、         研究反演算法（2.反演程序又和正演程序有什么关系？）

3、         实现反演算法并验证。（上面1.2两条是本人迷惘之处）

正演程序： dimension r(28),h(28),rr(50),ff(50),af(50)

complex z(28)

character*12 file1,file2,file3

pi=3.1415926

cc=0.0000001

c=8.0*pi**2*cc

write(*,*)\'input m=?  n=?  \'

read(*,*) m,n

do 5 i=1,n

write(*,*) \'number layar=\',i

write(*,*) \'input depth resistance\'

read(*,*) h(i),r(i)

5  continue

c write(*,*) \'please input freqence point file name\'

c read(*,\'(a)\')

file1=\'fre.dat\'

write(*,*)\'please input model file name\'

read(*,\'(a)\')file2

write(*,*)\'please input fc result file name\'

read(*,\'(a)\')file3

open(4,file=file2)

write(4,*)(h(i),r(i),i=1,n)

close(4)

open(8,file=file1,STATUS=\'OLD\')

read(8,*)(ff(i),i=1,m)

close(8)

do 15 i=1,m

f=ff(i)

call fc(r,h,z,n,f)

re=real(z(1))

aim=aimag(z(1))

af(i)=atan(aim/re)*180/pi

rr(i)=(re**2+aim**2)/(c*f)

15 continue

open(3,file=file3)

write(3,*) (1/ff(k),af(k),k=1,m)

close(3)

stop

end

subroutine fc(r,h,z,n,f)

dimension h(28),r(28)

complex z(28),a,b,d

pi=3.1415926

cc=0.0000001

c=8.0*pi**2*cc

z(n)=sqrt(cmplx(0.0,c*r(n)*f))

do 10 j=1,n-1

j1=n-j

a=sqrt(cmplx(0.0,c*f*r(j1)))

b=sqrt(cmplx(0.0,c*f/r(j1)))

d=exp(-2.0*b*h(j1))

z(j1)=a*(a*(1-d)+z(j1+1)*(1+d))/(a*(1+d)+z(j1+1)*(1-d))

10 continue

return

end



反演程序： dimension r(15),h(15),rr(50),ac(50,29),q(29),u(50,29),v(29,29),

     * e(29),ee(50),xz(29),ee1(50),r1(15),h1(15),ff(50),ft(29)

complex z(15),zz(15),zh(15),zr(15),z1r(15),z1h(15)

pi=3.1415926

cc=0.0000001

c=8.0*pi**2*cc

write(*,*)\'************************************************\'

write(*,*)\'**                                            **\'

write(*,*)\'**     m=?     n=?       lmax=?               **\'

write(*,*)\'**                                            **\'

write(*,*)\'************************************************\'

read(*,*) m,n,lmax

n1=2*n-1

open(8,file=\'fre.dat\',status=\'old\')

read(8,*) (ff(i),i=1,m)

        do 2 i=1,m

2 ff(i)=1.0/ff(i)

write(*,*) \'m=\',m, \'n=\',n, \'lmax=\',lmax

write(*,*)\'do you have field data if you have input 1\'

read(*,*) nn

if(nn.eq.1)  goto 123

do 5 i=1,n

write(*,*) \'number layar=\',i

write(*,*) \'input depth resistance\'

read(*,*) h(i),r(i)

5  continue

open(6,file=\'rr.dat\',status=\'new\')

do 15 i=1,m

f=ff(i)

call fc(r,h,z,n,f)

re=real(z(1))

aim=aimag(z(1))

rr(i)=(re**2+aim**2)/(c*f)

15 continue

write(6,*) (rr(i),i=1,m)

close(6)

123 continue

c------ initial model---------------------------------------

open(5,file=\'model\',status=\'old\')

read(5,*) (h(i),r(i),i=1,n)

close(5)

c-----------------------------------------------------------

c****** input field data ******************

open(6,file=\'rr.dat\',status=\'old\')

read(6,*) (rr(i),i=1,m)

close(6)

c*******************************************

sum=0.0

do 600 i=1,m

f=ff(i)

call fc(r,h,z,n,f)

s=1.0/(f*c)*((real(z(1)))**2+(aimag(z(1)))**2)

sum=sum+(s-rr(i))**2

ee(i)=rr(i)-s

600 continue

q0=sum

l=0

1 arf=1.e-10

  continue

l=l+1

write(*,*)\'l=\',l,  \'q0=\',q0

do 200 i=1,m

f=ff(i)

call fc(r,h,z,n,f)

call fc1(r,h,z,n,f,zz)

call fc2(r,h,z,zz,n,f,zh,zr)

call fc3(zz,zh,zr,n,f,z1h,z1r)

do 220 j=1,n

ac(i,j)=1.0/(c*f)*(z(1)*cmplx(real(z1r(j)),-aimag(z1r(j)))+

     * z1r(j)*cmplx(real(z(1)),-aimag(z(1))))

220 continue

do 230 j=1,n-1

ac(i,j+n)=1.0/(c*f)*(z(1)*cmplx(real(z1h(j)),-aimag(z1h(j)))+

     * z1h(j)*cmplx(real(z(1)),-aimag(z(1))))

230 continue

200 continue

call svd(m,n1,ac,q,u,v,e)

500  do 50 j=1,n1

yy=0.0

do 30 k=1,n1

if(q(k).le.td) goto 30

buk=0.0

do 40 i=1,m

buk=buk+ee(i)*u(i,k)

40 continue

qk=q(k)/(q(k)*q(k)+arf/(q(k)*q(k)))

yy=yy+qk*buk*v(j,k)

30 continue

xz(j)=yy

50 continue

do 650 i=1,n

650 r1(i)=r(i)+xz(i)

do 660 i=1,n-1

660 h1(i)=h(i)+xz(n+i)

sum1=0.0

do 680 i=1,m

f=ff(i)

call fc(r1,h1,z,n,f)

ss=1.0/(f*c)*((real(z(1)))**2+(aimag(z(1)))**2)

sum1=sum1+(rr(i)-ss)**2

ee1(i)=rr(i)-ss

680  continue

if(sum1.gt.q0) then

arf=arf*10.0

goto 500

else

do 690 j=1,n

r(j)=r1(j)

h(j)=h1(j)

690 continue

do 695 i=1,m

695 ee(i)=ee1(i)

q0=sum1

if(l.gt.lmax) goto 105

endif

c arf=arf/10.0

goto 1

105 write(*,*) \'iter gt lmax\'

open(7,file=\'result\',status=\'new\')

write(7,*) (h(i),r(i),i=1,n)

close(7)

stop

end

subroutine fc(r,h,z,n,f)

dimension h(15),r(15)

complex z(15),a,b,d

pi=3.1415926

cc=0.0000001

c=8.0*pi**2*cc

z(n)=sqrt(cmplx(0.0,c*r(n)*f))

do 10 j=1,n-1

j1=n-j

a=sqrt(cmplx(0.0,c*f*r(j1)))

b=sqrt(cmplx(0.0,c*f/r(j1)))

d=exp(-2.0*b*h(j1))

z(j1)=a*(a*(1-d)+z(j1+1)*(1+d))/(a*(1+d)+z(j1+1)*(1-d))

10 continue

return

end

subroutine fc1(r,h,z,n,f,zz)

complex zz(15),z(15),a,b,d,e

dimension r(15),h(15)

pi=3.1415926

cc=0.0000001

c=8.0*pi**2*cc

do 10 i=1,n-1

a=sqrt(cmplx(0.0,c*f*r(i)))

b=sqrt(cmplx(0.0,c*f/r(i)))

d=exp(-2.0*b*h(i))

e=a*(1+d)+z(i+1)*(1-d)

10 zz(i)=4.0*a**2/e**2*d

return

end

subroutine fc2(r,h,z,zz,n,f,zh,zr)

complex zz(15),z(15),zh(15),zr(15),a,b

dimension r(15),h(15)

pi=3.1415926

cc=0.0000001

c=8.0*pi**2*cc

a=cmplx(0.0,f*c)

b=sqrt(cmplx(0.0,c*f*r(n)))

do 10 i=1,n-1

10 zh(i)=(a-z(i+1)**2/r(i))*zz(i)

do 20 i=1,n-1

zr(i)=1.0/(2.0*r(i))*(z(i)-z(i+1)*zz(i)-h(i)*zh(i))

  20 continue

zr(n)=1.0/2.0*a/b

return

end

subroutine fc3(zz,zh,zr,n,f,z1h,z1r)

complex zz(15),zh(15),zr(15),z1h(15),z1r(15),x

do 10 i=1,n-1

x=cmplx(1.0,0.0)

if(i.eq.1) then

z1h(i)=zh(1)

z1r(i)=zr(1)

else

do 20 j=1,i-1

20 x=x*zz(j)

z1h(i)=x*zh(i)

z1r(i)=x*zr(i)

endif

10 continue

z1r(n)=x*zz(n-1)*zr(n)

return

end

C       THIS IS SVD PROGRAM

C

        SUBROUTINE SVD(M,N,A,Q,U,V,E)

        DIMENSION A(50,29),Q(29),U(50,29),V(29,29),E(29)

        LOGICAL WITHU,WITHV

        EPS=1.E-14

        TOL=1.E-16

        WITHU=.TRUE.

        WITHV=.TRUE.

        DO 10 I=1,M

        DO 10 J=1,N

10      U(I,J)=A(I,J)

        G=0.

        X=0.

C       HOUSHOLDER\'S REDUCTION TO BIDIAGNAL FORM

        DO 200 I=1,N

        E(I)=G

        S=0.

        L=I+1

        DO 20 J=I,M

20      S=S+U(J,I)**2

        IF(S.GE.TOL)GOTO 30

        G=0.

        GOTO 90

30      F=U(I,I)

        IF(F.LT.0.)GOTO 40

        G=-SQRT(S)

        GOTO 50

40      G=SQRT(S)

50      H=F*G-S

        U(I,I)=F-G

        IF(L.GT.N)GOTO 90

        DO 80 J=L,N

        S=0.

        DO 60 K=I,M

60      S=S+U(K,I)*U(K,J)

        F=S/H

        DO 70 K=I,M

70      U(K,J)=U(K,J)+F*U(K,I)

80      CONTINUE

90      Q(I)=G

        S=0.

        IF(L.GT.N)GOTO 102

        DO 100 J=L,N

100     S=S+U(I,J)**2

102     IF(S.GE.TOL)GOTO 110

        G=0.0

        GOTO 180

110     F=U(I,I+1)

        IF(F.LT.0.)GOTO 120

        G=-SQRT(S)

        GOTO 130

120     G=SQRT(S)

130     H=F*G-S

        U(I,I+1)=F-G

        DO 140 J=L,N

140     E(J)=U(I,J)/H

        DO 170 J=L,M

        S=0.

        DO 150 K=L,N

150     S=S+U(J,K)*U(I,K)

        DO 160 K=L,N

160     U(J,K)=U(J,K)+S*E(K)

170     CONTINUE

180     Y=ABS(Q(I))+ABS(E(I))

        IF(Y.LE.X)GOTO 200

        X=Y

200     CONTINUE

C       ACCUMULATION OF RIGHT-HAND TRANSFORMATION

        IF(WITHV)GO TO 210

        GOTO 290

210     DO 280 II=1,N

        I=N+1-II

        IF(G.EQ.0.)GOTO 260

        H=U(I,I+1)*G

        DO 220 J=L,N

220     V(J,I)=U(I,J)/H

        DO 250 J=L,N

        S=0.

        DO 230 K=L,N

230     S=S+U(I,K)*V(K,J)

        DO 240 K=L,N

240     V(K,J)=V(K,J)+S*V(K,I)

250     CONTINUE

260     IF(L.GT.N)GOTO 272

        DO 270 J=L,N

        V(I,J)=0.

270     V(J,I)=0.

272     V(I,I)=1.

        G=E(I)

        L=I

280     CONTINUE

C       ACCUMULATION OF LEFT-HANDM TRANSFORMATION

290     IF(WITHU)GOTO 300

        GOTO 400

300     DO 390 II=1,N

        I=N+1-II

        L=I+1

        G=Q(I)

        IF(L.GT.N)GOTO 311

        DO 310 J=L,N

310     U(I,J)=0.

311     IF(G.EQ.0.)GOTO 360

        H=U(I,I)*G

        IF(L.GT.N)GOTO 341

        DO 340 J=L,N

        S=0.

        DO 320 K=L,M

320     S=S+U(K,I)*U(K,J)

        F=S/H

        DO 330 K=I,M

330     U(K,J)=U(K,J)+F*U(K,I)

340     CONTINUE

341     DO 350 J=I,M

350     U(J,I)=U(J,I)/G

        GOTO 380

360     DO 370 J=I,M

370     U(J,I)=0.

380     U(I,I)=U(I,I)+1.

390     CONTINUE

C       DIAGONALIZATION OF THE BIDIAGONAL FORM

400     EPS=EPS*X

        DO 580 KK=1,N

        K=N+1-KK

410     DO 420 LL=1,K

        L=K+1-LL

        IF(ABS(E(L)).LE.EPS)GOTO 470

        IF(ABS(Q(L-1)).LE.EPS)GOTO 430

420     CONTINUE

C       CANCELLATION OF E(L) IF L>1

430     C=0.

        S=1.

        L1=L-1

        DO 460 I=L,K

        F=E(I)*S

        E(I)=C*E(I)

        IF(ABS(F).LE.EPS)GOTO 470

        G=Q(I)

        H=SQRT(F*F+G*G)

        Q(I)=H

        C=G/H

        S=-F/H

        IF(WITHU) GOTO 440

        GOTO 460

440     DO 450 J=1,M

        Y=U(J,L1)

        Z=U(J,I)

        U(J,L1)=Y*C+Z*S

450     U(J,I)=-Y*S+Z*C

460     CONTINUE

470     Z=Q(K)

        IF(L.EQ.K)GOTO 550

C       SHIFT FORM BOTTOM IXI MINOR

        X=Q(L)

        Y=Q(K-1)

        G=E(K-1)

        H=E(K)

        F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.*H*Y)

        G=SQRT(F*F+1.)

        IF(F.LT.0.)GOTO 480

        F=((X-Z)*(X+Z)+H*(Y/(F+G)-H))/X

        GOTO 490

480     F=((X-Z)*(X+Z)+H*(Y/(F-G)-H))/X

C       NEXT QR TRANSFORMATION

490     S=1.

        C=1.

        LO=L+1

        DO 540 I=LO,K

        G=E(I)

        Y=Q(I)

        H=S*G

        G=C*G

        E(I-1)=SQRT(F*F+H*H)

        Z=E(I-1)

        C=F/Z

        S=H/Z

        F=X*C+G*S

        G=-X*S+G*C

        H=Y*S

        Y=Y*C

        IF(WITHV)GOTO 495

        GOTO 510

495     DO 500 J=1,N

        X=V(J,I-1)

        Z=V(J,I)

        V(J,I-1)=X*C+Z*S

500     V(J,I)=-X*S+Z*C

510     Q(I-1)=SQRT(F*F+H*H)

        Z=Q(I-1)

        C=F/Z

        S=H/Z

        F=C*G+S*Y

        X=-S*G+C*Y

        IF(WITHU)GOTO 520

        GOTO 540

520     DO 530 J=1,M

        Y=U(J,I-1)

        Z=U(J,I)

        U(J,I-1)=Y*C+Z*S

530     U(J,I)=-Y*S+Z*C

540     CONTINUE

        E(L)=0.

        E(K)=F

        Q(K)=X

        GOTO 410

550     IF(Z.GE.0.)GOTO 580

C       Q(K) IS MODE NON-NEGATIVE

        Q(K)=-Z

        IF(WITHV) GOTO 560

        GO TO 580

560     DO 570 J=1,N

570     V(J,K)=-V(J,K)

580     CONTINUE

c       WRITE(*,581)

581     FORMAT(1X,\'   SINGULE----VALUE\')

c       WRITE(*,590)(Q(II),II=1,N)

590     FORMAT(1X,8(2X,E13.7))

        RETURN

        END

xiexie !