虚边界元方法：圆域位势问题求解程序

C=============================虚边界  圆域 线性 求内点及边界点的U值===============================
C====================================MAIN============================================      

DIMENSION XP(500),YP(500),CX(50),CY(50),C(500),SC(500)
DIMENSION KODE(500),FI(500),SOL(50),F(500)
DIMENSION H(500,500),G(500,500),DFI(500)
OPEN(11,FILE='E:\XUBIANJIE\XU\.TXT',STATUS='OLD')
OPEN(16,FILE='E:\XUBIANJIE\XU\.OUT',STATUS='OLD')
CALL INPUT(ARF,R,D,NE,NIT,SC)
CALL PREPDATA(NE,NIT,R,D,C,SC,XP,YP,CX,CY,FI,KODE,ARF)
CALL CHG(NE,C,R,D,XP,YP,H,G,KODE,ARF)
CALL SEQ(NE,H,G,DFI,FI,F,KODE)
CALL INTER(R,D,C,NE,NIT,H,G,F,CX,CY,SOL,ARF)
CALL OUTPUT(NE,NIT,C,R,XP,YP,FI,DFI,CX,CY,SOL)
STOP
END
C==========================INPUT====================================================
      SUBROUTINE INPUT(ARF,R,D,NE,NIT,SC)
DIMENSION SC(50)
READ(11,*) ARF,R,D,NE,NIT
READ(11,*)(SC(I),I=1,NIT)
      RETURN
END
C=========================PREPDATA================================================
      SUBROUTINE PREPDATA(NE,NIT,R,D,C,SC,XP,YP,CX,CY,FI,KODE,ARF)
 DIMENSION XP(2*NE),YP(2*NE),C(NE+1),
     +FI(500),KODE(500),CX(NIT),CY(NIT),SC(NIT)
      TWOPAI=ATAN(1.0)*8.0
C(NE+1)=TWOPAI
DO 10 I=1,NE

C(I)=TWOPAI*(I-1)/NE
10 CONTINUE
DO 20 I=1,NE
XP(2*I-1)=R*DCOS((1+ARF)*C(I)/2.0+(1-ARF)*C(I+1)/2.0)
      XP(2*I)=R*DCOS((1-ARF)*C(I)/2.0+(1+ARF)*C(I+1)/2.0)
YP(2*I-1)=R*DSIN((1+ARF)*C(I)/2.0+(1-ARF)*C(I+1)/2.0)
      YP(2*I)=R*DSIN((1-ARF)*C(I)/2.0+(1+ARF)*C(I+1)/2.0)
20 CONTINUE
      DO 50 I=1,NE
  KODE(I)=0.0
  FI(I)=YP(I)
50 CONTINUE
      DO 60 I=NE+1,2*NE
  KODE(I)=0.0
  FI(I)=YP(I)
60 CONTINUE
      DO 70 I=1,NIT
  SC(I)=SC(I)/360*TWOPAI
  CX(I)=(R/2)*DCOS(SC(I))
  CY(I)=(R/2)*DSIN(SC(I))
70 CONTINUE
RETURN
END
C==========================================CHG===============================
      SUBROUTINE CHG(NE,C,R,D,XP,YP,H,G,KODE,ARF)
DIMENSION C(NE+1),XP(2*NE),YP(2*NE),H(500,500),G(500,500),KODE(500)
      DO 5 I=1,2*NE
DO 5 J=1,2*NE
      H(I,J)=0.0
G(I,J)=0.0