C /SUBA/ SUB1AHL SUBROUTINE SUBROUTINE GRAP7G(IA,LB,LL,NX,NXX,NXP,X0,Y0,XL,YL, & U,VMIN,VMAX,KP) DIMENSION U(1) C XX=999.0 XM=0.075 YM=0.075 C NX1=NX+1 NX2=NX+2 N1=NXX+2 C INX=NXX/NX INX=1 L1=LB+LL-1 DO 8 J=LB,L1 DO 8 I=1,NX2 I1=I+N1*(J-1) I2=I*INX+N1*(J-1) U(I1)=U(I2) 8 CONTINUE DX=XL/FLOAT(NX1) I1=1+N1*(LB-1) UMAX=U(I1)+1.0E-6 UMIN=U(I1)-1.0E-6 DO 10 J=LB,L1 DO 10 I=1,NX2 I1=I+N1*(J-1) UMAX=AMAX1(UMAX,U(I1)) UMIN=AMIN1(UMIN,U(I1)) 10 CONTINUE IF(KP.EQ.0) GO TO 16 IF(KP.EQ.1) GO TO 12 IF(KP.EQ.-1) GO TO 14 IF(KP.EQ.5) GO TO 15 UMIN=VMIN UMAX=VMAX GO TO 16 12 UMAX=VMAX GO TO 16 14 UMIN=VMIN GO TO 16 15 UMIN=-UMIN UMAX=AMAX1(UMAX,UMIN) UMIN=-UMAX 16 CONTINUE C DO 18 J=LB,L1 DO 18 I=1,NX2 I1=I+N1*(J-1) IF(U(I1).GT.UMAX) U(I1)=UMAX IF(U(I1).LT.UMIN) U(I1)=UMIN 18 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,1,IA,UMIN,UMAX) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) C X=X0+0.5*XL*FLOAT(NX1-2*NXP)/FLOAT(NX1) Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) call linee C DO 20 JJ=LB,L1 I=JJ+IA-LB-1 C I=3-MOD(JJ-LB,3) I=3-MOD(I,3) CALL NEWPE1(I) X=X0 I1=1+N1*(JJ-1) Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,3) DO 30 I=1,NX I1=I+1+N1*(JJ-1) X=X+DX Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,2) 30 CONTINUE I1=NX2+N1*(JJ-1) X=X+DX Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,2) call linee 20 CONTINUE RETURN END SUBROUTINE GRAP4M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX, & KP,LANK,V0) C DIMENSION U(1) XX=999.0 XM=0.075 YM=0.075 NXY=NX*NY NX1=NX+1 NY1=NY+1 NX2=NX+2 NY2=NY+2 N1=NX2 N2=N1*NY2 C DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 10 J=1,NY2 DO 10 I=1,NX2 Y=0.5*HY*FLOAT(2*J-NY-3) X=0.5*HX*FLOAT(2*I-NX-3+2*NXP) R=X*X+Y*Y IF(R.LT.AR2) GO TO 10 I1=I+NX2*(J-1) X1=U(I1) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 10 CONTINUE IF(KP.EQ.0) GO TO 16 IF(KP.EQ.1) GO TO 12 IF(KP.EQ.-1) GO TO 14 IF(KP.EQ.5) GO TO 15 UMIN=VMIN UMAX=VMAX GO TO 16 12 UMAX=VMAX GO TO 16 14 UMIN=VMIN GO TO 16 15 CONTINUE X=ABS(UMAX-V0) Y=ABS(V0-UMIN) X=AMAX1(X,Y) UMAX=V0+X UMIN=V0-X 16 CONTINUE C DO 18 I=1,N2 IF(U(I).GT.UMAX) U(I)=UMAX IF(U(I).LT.UMIN) U(I)=UMIN 18 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,KP,UMIN,UMAX) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) Y=Y0+0.5*YL CALL PLOT(X,Y,3) X=X0+XL CALL PLOT(X,Y,2) X=X0+XL*(0.5-FLOAT(NXP)/FLOAT(NX1)) Y=Y0 CALL PLOT(X,Y,3) Y=Y+YL CALL PLOT(X,Y,2) C call linee C ULD=FLOAT(LANK+1)/(UMAX-UMIN) EP1=1.0E-8 IPE1=0 C DO 110 J=1,NY1 CALL NEWPE1(IPE1) DO 100 I=1,NX1 C I1=I+N1*(J-1) I2=I1+1 I3=I1+1+N1 I4=I1+N1 C U1=(U(I1)-UMIN)*ULD U2=(U(I2)-UMIN)*ULD U3=(U(I3)-UMIN)*ULD U4=(U(I4)-UMIN)*ULD C K1=IFIX(U1) K2=IFIX(U2) K3=IFIX(U3) K4=IFIX(U4) C J1=IABS(K2-K1) J2=IABS(K3-K2) J3=IABS(K4-K3) J4=IABS(K1-K4) C IF(J1.EQ.0) GO TO 28 DO 20 I1=1,J1 U0=FLOAT(I1)+FLOAT(MIN0(K1,K2)) IF(ABS(U2-U1).LT.EP1) GO TO 20 X1=DX*((U0-U1)/(U2-U1)+FLOAT(I-1)) Y1=DY*FLOAT(J-1) X1=X0+X1 Y1=Y0+Y1 C IF((U3-U0)*(U2-U0).GT.0.0) GO TO 22 IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 22 IF(ABS(U3-U2).LT.EP1) GO TO 22 X2=DX*FLOAT(I) Y2=DY*((U0-U2)/(U3-U2)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 22 CONTINUE C IF((U4-U0)*(U3-U0).GT.0.0) GO TO 24 IF((U1-U0)*(U3-U0).GT.0.0) GO TO 24 IF((U2-U0)*(U4-U0).GT.0.0) GO TO 24 IF(ABS(U4-U3).LT.EP1) GO TO 24 X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y2=DY*FLOAT(J) X2=X0+X2 Y2=Y0+Y2 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 24 CONTINUE C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 26 IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 26 IF(ABS(U1-U4).LT.EP1) GO TO 26 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 26 CONTINUE 20 CONTINUE 28 CONTINUE C IF(J2.EQ.0) GO TO 38 DO 30 I1=1,J2 U0=FLOAT(I1)+FLOAT(MIN0(K2,K3)) IF(ABS(U3-U2).LT.EP1) GO TO 30 X1=DX*FLOAT(I) Y1=DY*((U0-U2)/(U3-U2)+FLOAT(J-1)) X1=X0+X1 Y1=Y0+Y1 C IF((U4-U0)*(U3-U0).GT.0.0) GO TO 32 IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 32 IF(ABS(U4-U3).LT.EP1) GO TO 32 X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y2=DY*FLOAT(J) X2=X0+X2 Y2=Y0+Y2 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 32 CONTINUE C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 34 IF((U1-U0)*(U3-U0).GT.0.0) GO TO 34 IF((U2-U0)*(U4-U0).GT.0.0) GO TO 34 IF(ABS(U1-U4).LT.EP1) GO TO 34 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 34 CONTINUE 30 CONTINUE 38 CONTINUE C IF(J3.EQ.0) GO TO 48 DO 40 I1=1,J3 U0=FLOAT(I1)+FLOAT(MIN0(K3,K4)) IF(ABS(U4-U3).LT.EP1) GO TO 40 X1=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y1=DY*FLOAT(J) X1=X0+X1 Y1=Y0+Y1 C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 42 IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 42 IF(ABS(U1-U4).LT.EP1) GO TO 42 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 42 CONTINUE 40 CONTINUE 48 CONTINUE C 100 CONTINUE call linee 110 CONTINUE RETURN END SUBROUTINE GRAP4N(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX, & KP,LANK,V0,TH0,DTH,IPEN,AFAC) C DIMENSION U(1) PI=3.1415926 XX=999.0 XM=0.075 YM=0.075 NXY=NX*NY NX1=NX+1 NY1=NY+1 NX2=NX+2 NY2=NY+2 N1=NX2 N2=N1*NY2 IHI=3 ILO=2 CALL NEWPE1(IPEN) C DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 10 J=1,NY2 DO 10 I=1,NX2 Y=0.5*HY*FLOAT(2*J-NY-3) X=0.5*HX*FLOAT(2*I-NX-3+2*NXP) R=X*X+Y*Y IF(R.LT.AR2) GO TO 10 I1=I+NX2*(J-1) X1=U(I1) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 10 CONTINUE IF(KP.EQ.0) GO TO 16 IF(KP.EQ.1) GO TO 12 IF(KP.EQ.-1) GO TO 14 IF(KP.EQ.5) GO TO 15 UMIN=VMIN UMAX=VMAX GO TO 16 12 UMAX=VMAX GO TO 16 14 UMIN=VMIN GO TO 16 15 CONTINUE X=ABS(UMAX-V0) Y=ABS(V0-UMIN) X=AMAX1(X,Y) UMAX=V0+X UMIN=V0-X 16 CONTINUE C DO 18 I=1,N2 IF(U(I).GT.UMAX) U(I)=UMAX IF(U(I).LT.UMIN) U(I)=UMIN 18 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,KP,UMIN,UMAX) C CALL NEWPE1(ILO) CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) Y=Y0+0.5*YL CALL PLOT(X,Y,3) X=X0+XL CALL PLOT(X,Y,2) X=X0+0.5*XL Y=Y0 CALL PLOT(X,Y,3) Y=Y0+YL CALL PLOT(X,Y,2) C X1=COS(PI*TH0/180.0) X=X0+0.5*XL Y=Y0+0.5*YL DO 50 I=1,20 TH=TH0+DTH*FLOAT(I-1) R=PI*TH/180.0 R=0.5*XL*COS(R)/X1 IF(R.LE.0.0) GO TO 52 CALL GCIRC1(X,Y,R,90) 50 CONTINUE 52 CONTINUE call linee C ULD=FLOAT(LANK+1)/(UMAX-UMIN) U00=(V0-UMIN)*ULD EP1=1.0E-8 XO0=X0+0.5*XL YO0=Y0+0.5*YL RO0=(0.55*XL)**2 IPE1=0 C DO 110 J=1,NY1 c CALL NEWPE1(IPE1) DO 100 I=1,NX1 C I1=I+N1*(J-1) I2=I1+1 I3=I1+1+N1 I4=I1+N1 C U1=(U(I1)-UMIN)*ULD U2=(U(I2)-UMIN)*ULD U3=(U(I3)-UMIN)*ULD U4=(U(I4)-UMIN)*ULD X1=X0+DX*FLOAT(I-1) X2=X1+DX Y1=Y0+DY*FLOAT(J-1) Y2=Y1+DY U0=0.25*(U1+U2+U3+U4) X=U00+AFAC*(UMAX-UMIN)*ULD IF(U0.LT.X) GO TO 60 c CALL NEWPE1(1) c CALL PLOT(X1,Y1,3) c CALL PLOT(X2,Y2,2) c CALL PLOT(X1,Y2,3) c CALL PLOT(X2,Y1,2) 60 CONTINUE C K1=IFIX(U1) K2=IFIX(U2) K3=IFIX(U3) K4=IFIX(U4) C J1=IABS(K2-K1) J2=IABS(K3-K2) J3=IABS(K4-K3) J4=IABS(K1-K4) C IF(J1.EQ.0) GO TO 28 DO 20 I1=1,J1 U0=FLOAT(I1)+FLOAT(MIN0(K1,K2)) IF(U0.GE.U00) CALL NEWPE1(IHI) IF(U0.LT.U00) CALL NEWPE1(ILO) IF(ABS(U2-U1).LT.EP1) GO TO 20 X1=DX*((U0-U1)/(U2-U1)+FLOAT(I-1)) Y1=DY*FLOAT(J-1) X1=X0+X1 Y1=Y0+Y1 C IF((U3-U0)*(U2-U0).GT.0.0) GO TO 22 IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 22 IF(ABS(U3-U2).LT.EP1) GO TO 22 X2=DX*FLOAT(I) Y2=DY*((U0-U2)/(U3-U2)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 22 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 22 CONTINUE C IF((U4-U0)*(U3-U0).GT.0.0) GO TO 24 IF((U1-U0)*(U3-U0).GT.0.0) GO TO 24 IF((U2-U0)*(U4-U0).GT.0.0) GO TO 24 IF(ABS(U4-U3).LT.EP1) GO TO 24 X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y2=DY*FLOAT(J) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 24 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 24 CONTINUE C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 26 IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 26 IF(ABS(U1-U4).LT.EP1) GO TO 26 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 26 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 26 CONTINUE 20 CONTINUE 28 CONTINUE C IF(J2.EQ.0) GO TO 38 DO 30 I1=1,J2 IF(U0.GE.U00) CALL NEWPE1(IHI) IF(U0.LT.U00) CALL NEWPE1(ILO) U0=FLOAT(I1)+FLOAT(MIN0(K2,K3)) IF(ABS(U3-U2).LT.EP1) GO TO 30 X1=DX*FLOAT(I) Y1=DY*((U0-U2)/(U3-U2)+FLOAT(J-1)) X1=X0+X1 Y1=Y0+Y1 C IF((U4-U0)*(U3-U0).GT.0.0) GO TO 32 IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 32 IF(ABS(U4-U3).LT.EP1) GO TO 32 X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y2=DY*FLOAT(J) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 32 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 32 CONTINUE C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 34 IF((U1-U0)*(U3-U0).GT.0.0) GO TO 34 IF((U2-U0)*(U4-U0).GT.0.0) GO TO 34 IF(ABS(U1-U4).LT.EP1) GO TO 34 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 34 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 34 CONTINUE 30 CONTINUE 38 CONTINUE C IF(J3.EQ.0) GO TO 48 DO 40 I1=1,J3 U0=FLOAT(I1)+FLOAT(MIN0(K3,K4)) IF(U0.GE.U00) CALL NEWPE1(IHI) IF(U0.LT.U00) CALL NEWPE1(ILO) IF(ABS(U4-U3).LT.EP1) GO TO 40 X1=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y1=DY*FLOAT(J) X1=X0+X1 Y1=Y0+Y1 C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 42 IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 42 IF(ABS(U1-U4).LT.EP1) GO TO 42 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 42 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 42 CONTINUE 40 CONTINUE 48 CONTINUE C 100 CONTINUE call linee 110 continue CALL NEWPE1(IPEN) RETURN END SUBROUTINE GRAP3A(NR,NZ,MB,NX,NY,NXP,X0,Y0,XL,YL, & AR2,U,VMIN,VMAX,URMI0,KP,IAR,LAS1) PARAMETER (XM=3.4, YM=-0.6) DIMENSION U(1) C CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8 C NR1=NR+1 NR2=NR+2 NZ1=NZ+1 NZ2=NZ+2 N1=NR2 N2=NR2*NZ2 HX=FLOAT(NR+1)/FLOAT(NX+1) HY=FLOAT(NZ+1)/FLOAT(NY+1) RMIN=0.0 RMAX=FLOAT(NR1) ZMIN=0.0 ZMAX=FLOAT(NZ+1) C NXY=NX*NY NL=NXY*2 AD=3.0*HX/FLOAT(LAS1) AD=0.3*HX AD=0.1*HX DX=XL/FLOAT(NR+1) DY=YL/FLOAT(NZ+1) AR=DX*HX*0.20 CR=0.40 C UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 8 J=1,NZ2 DO 8 I=1,NR2 Y=0.5*DY*FLOAT(2*J-NZ2-1) X=0.5*DX*FLOAT(2*I-NR2-1+2*NXP) R=X*X+Y*Y C IF(R.LT.AR2) GO TO 8 I1=I+N1*(J-1) I2=I1+N2 X1=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 8 CONTINUE IF(KP.EQ.0) GO TO 126 IF(KP.EQ.1) GO TO 122 IF(KP.EQ.-1) GO TO 124 UMIN=VMIN UMAX=VMAX GO TO 126 122 UMAX=VMAX GO TO 126 124 UMIN=VMIN 126 CONTINUE C UMIN=-UMIN UMAX=AMAX1(UMAX,UMIN) UMIN=-UMAX URMI=URMI0*UMAX C DO 18 I=1,N2 I1=I+N2*(MB-1) I2=I1+N2 R=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) IF(R.GT.UMAX) U(I1)=U(I1)*UMAX/R IF(R.GT.UMAX) U(I2)=U(I2)*UMAX/R 18 CONTINUE C C ENCODE(M3,100) UMIN C ENCODE(M4,100) UMAX 100 FORMAT(1PE8.1) X1=X0+XM Y1=Y0+YM C CALL SYMBOL(X0,Y1,0.25,M1,0.0,4) C CALL SYMBOL(999.,Y1,0.25,M3,0.0,8) C CALL SYMBOL(X1,Y1,0.25,M2,0.0,4) C CALL SYMBOL(999.,Y1,0.25,M4,0.0,8) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) X=X0+0.5*XL*FLOAT(NR1-2*NXP)/FLOAT(NR1) Y=Y0 CALL PLOT(X,Y,3) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 Y=Y0+0.5*YL CALL PLOT(X,Y,3) X=X0+XL CALL PLOT(X,Y,2) call linee IPE1=0 C NY1=NY/2+1 DO 10 J=1,NY CALL NEWPE1(IPE1) DO 12 I=1,NX R1=FLOAT(I)*HX Z1=FLOAT(J)*HY X1=X0+DX*R1 Y1=Y0+DY*Z1 X2=X1 Y2=Y1 RR=0 C DO 20 III=1,LAS1 X1=X2 Y1=Y2 IR1=IFIX(R1)+1 IZ1=IFIX(Z1)+1 IF(IR1.LT.1.OR.IR1.GT.NR1) GO TO 20 IF(IZ1.LT.1.OR.IZ1.GT.NZ1) GO TO 20 I1=IR1+NR2*(IZ1-1) I2=I1+N2 R=R1-IFIX(R1) Z=Z1-IFIX(Z1) UX=R*Z*U(I1+N1+1)+(1.0-R)*Z*U(I1+N1) & +R*(1.0-Z)*U(I1+1)+(1.0-R)*(1.0-Z)*U(I1) UY=R*Z*U(I2+N1+1)+(1.0-Z)*R*U(I2+N1) & +R*(1.0-Z)*U(I2+1)+(1.0-R)*(1.0-Z)*U(I2) R=SQRT(UX*UX+UY*UY) IF(R.LT.0.01*URMI) GO TO 12 IF(R.GT.UMAX) UX=UX*UMAX/R IF(R.GT.UMAX) UY=UY*UMAX/R UX=AD*UX/UMAX UY=AD*UY/UMAX IF(R.LT.UMAX) RR=RR+R/UMAX IF(R.GE.UMAX) RR=RR+1.0 R2=R1+UX Z2=Z1+UY X2=X0+DX*R2 Y2=Y0+DY*Z2 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) R1=R2 Z1=Z2 IF(R.LT.URMI) GO TO 22 IF(R1.LE.RMIN.OR.R1.GE.RMAX) GO TO 22 IF(Z1.LE.ZMIN.OR.Z1.GE.ZMAX) GO TO 22 20 CONTINUE 22 CONTINUE XX1=RR/FLOAT(LAS1) AR1=AR IF(IAR.EQ.1) AR1=AR*XX1*1.5 CALL GAROHD(X1,Y1,X2,Y2,AR1,CR) 12 CONTINUE call linee 10 CONTINUE RETURN END SUBROUTINE BOUND3(NX,NY,NZ,NXP,MDFU,NB,BIS,RO01,PR01,VSW, & X1,HX,HY,HZ,ARU,AR2,F,P) DIMENSION F(1),P(1) NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 NZ1=NZ+1 NZ2=NZ+2 N1=NX2 N2=N1*NY2 N3=N2*NZ2 C NX11=NX1-2*NXP DO 12 K=1,NZ2 DO 12 J=1,NY2 DO 12 I=1,NX2 X=FLOAT(2*I-2-NX11)/FLOAT(NX11) X=(X1+1.0)*X+X1 IF(X.GT.0.0) X=0.0 X=X*X I1=I+N1*(J-1)+N2*(K-1) I2=I1+N3 I3=I2+N3 I4=I3+N3 I5=I4+N3 I6=I5+N3 I7=I6+N3 I8=I7+N3 F(I1)=(1.0-X)*F(I1)+RO01*X F(I2)=(1.0-X)*F(I2)+VSW*X F(I3)=(1.0-X)*F(I3) F(I4)=(1.0-X)*F(I4) F(I5)=(1.0-X)*F(I5)+PR01*X F(I6)=(1.0-X)*F(I6) F(I7)=(1.0-X)*F(I7) F(I8)=(1.0-X)*F(I8)+BIS*X 12 CONTINUE C DO 410 II1=1,MDFU C BOUNDARY CONDITION AT NZ=1 AND NZ=NZ2 C BOUNDARY CONDITION AT NY=1 AND NY=NY2 DO 110 M=1,NB DO 111 J=2,NY1 DO 111 I=2,NX1 I1=I+N1*(J-1)+N3*(M-1) I2=I1+N2*NZ1 I3=N2*2 IF(I.EQ.2) F(I2)=F(I2-N2-1) IF(I.GT.2) F(I2)=(4.0*F(I2-N2-1)-F(I2-I3-2))/3.0 F(I1)=F(I1+N2) IF(M.EQ.4) F(I1)=-F(I1) IF(M.EQ.6) F(I1)=-F(I1) IF(M.EQ.7) F(I1)=-F(I1) 111 CONTINUE DO 112 K=1,NZ2 DO 112 I=2,NX1 I1=I+N2*(K-1)+N3*(M-1) I2=I1+N1*NY1 I3=N1*2 IF(I.EQ.2) F(I2)=F(I2-N1-1) IF(I.GT.2) F(I2)=(4.0*F(I2-N1-1)-F(I2-I3-2))/3.0 F(I1)=F(I1+N1) IF(M.EQ.3) F(I1)=-F(I1) IF(M.EQ.7) F(I1)=-F(I1) 112 CONTINUE C BOUNDARY CONDITION AT NX=NX2 DO 113 K=1,NZ2 DO 113 J=1,NY2 I1=1+N1*(J-1)+N2*(K-1)+N3*(M-1) I2=I1+NX1 F(I2)=(4.0*F(I2-1)-F(I2-2))/3.0 113 CONTINUE 110 CONTINUE IF(II1.EQ.MDFU) GO TO 410 DO 414 M=1,NB DO 414 K=2,NZ1 L1=N2*(K-1)+N3*(M-1) DO 412 I=1,N2 I1=I+L1 P(I+2*N2)=F(I1+N2) P(I+N2)=F(I1) IF(K.EQ.2) P(I)=F(I1-N2) 412 CONTINUE DO 414 J=2,NY1 DO 414 I=2,NX1 J1=I+N1*(J-1)+L1 I1=I+N1*(J-1)+N2 F(J1)=(6.0*P(I1)+P(I1-1)+P(I1+1)+P(I1-N1)+P(I1+N1) & +P(I1-N2)+P(I1+N2))/12.0 P(I1-N2)=P(I1) 414 CONTINUE 410 CONTINUE C DO 14 K=1,NZ2 DO 14 J=1,NY2 DO 14 I=1,NX2 X=0.5*HX*FLOAT(2*I-NX1-2+2*NXP) Y=0.5*HY*FLOAT(2*J-3) Z=0.5*HZ*FLOAT(2*K-3) RO2=X*X+Y*Y+Z*Z I1=I+N1*(J-1)+N2*(K-1) I2=I1+N3 I3=I2+N3 I4=I3+N3 I5=I4+N3 I6=I5+N3 I7=I6+N3 I8=I7+N3 C X2=RO2/AR2-1.0 IF(X2.LT.0.0) X2=0.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) F(I2)=F(I2)*X2 F(I3)=F(I3)*X2 F(I4)=F(I4)*X2 14 CONTINUE C RETURN END SUBROUTINE AINTE1(IA,AA,F,P) DIMENSION IA(1),AA(1),F(1),P(1),Q(16) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO C AA 1 2 3 4 5 6 7 8 9 10 C TH0 ARU AR1 ARB XXL YYL ZZL B0 DL SIGN C 11 12 13 14 15 16 17 18 19 20 C DIRE BOUN GX0 GY0 GXL GYL GTH GXMI GXMA EP1 C P(1)=P(1) PI=3.1415926 NX=IA(1) NY=IA(2) MX=IA(3) MY=IA(4) MZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) TH0=AA(1)*PI/180.0 ARU=AA(2) AR1=AA(3) ARB=AA(4) XXL=AA(5) YYL=AA(6) ZZL=AA(7) B0=AA(8) GX0=AA(13) GY0=AA(14) GXL=AA(15) GYL=AA(16) GTH=AA(17)*PI/180.0 GXMI=AA(18) GXMA=AA(19) EP1=AA(20) C XX=999.0 XM=0.075 YM=0.075 X1=GX0 Y1=GY0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,MI,GXMI,GXMA) C XC1=COS(GTH) XS1=SIN(GTH) Y=GY0+0.5*GYL CALL PLOT(GX0,Y,3) X=GX0 Y=GY0+GYL CALL PLOT(X,Y,2) X=GX0+GXL CALL PLOT(X,Y,2) Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=GX0 CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1 Y=GY0+0.5*GYL*XS1+0.5*GYL CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1+GXL CALL PLOT(X,Y,2) X=GX0+GXL Y=GY0+0.5*GYL CALL PLOT(X,Y,2) C XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+GYL CALL PLOT(X,Y,3) Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=X+0.5*GYL*XC1 Y=Y+0.5*GYL*XS1 CALL PLOT(X,Y,2) C XL0=COS(TH0) XL01=XL0*1.0 XL=2.0*XL0 YL=XL NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 N1=NX2 N2=N1*NY2 AR2=AR1*AR1 DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) HX=XXL/FLOAT(MX+1) HY=YYL/FLOAT(MY+1) HZ=ZZL/FLOAT(MZ+1) XLMI=0.5*HX*FLOAT(2*NXP-MX-1) XLMA=0.5*HX*FLOAT(2*NXP+MX+1) DL0=HX*AA(9) DL=DL0 LAST=IFIX(XXL/DL) GYMI=-0.1*HY GYMA=AA(21) AMPY=0.5*(GXMA-GXMI)*GYL/(GYMA*GXL) GZMI=-0.1*HZ GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 C JB1=NY2/2+1 DO 100 J=JB1,NY2 DO 100 I=1,NX2 X=0.5*DX*(2*I-NX2-1) Y=0.5*DY*(2*J-NY2-1) R=SQRT(X*X+Y*Y) IF(R.GT.XL01) GO TO 100 Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=ARB*ARB-X1*X1-Y1*Y1 IF(Z1.LE.0.0) GO TO 100 Z1=SQRT(Z1) CALL QUANT1(X1,Y1,Z1,HX,HY,HZ,1,MX,MY,MZ,NXP,AA,F,Q) B1=SQRT(Q(M1)*Q(M1)+Q(M2)*Q(M2)+Q(M3)*Q(M3)) IF(B1.LT.1.0E-20) GO TO 100 AA(11)=1.0 X2=Q(M1)*X1+Q(M2)*Y1+Q(M3)*Z1 IF(X2.LT.0.0) AA(11)=-1.0 C X2=0.0 Y2=0.0 Z2=0.0 X3=0.0 Y3=0.0 Z3=0.0 GX1=GX0+XDD1+AMP*(X1+AMPY*Y1*XC1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) CALL PLOT(GX1,GY1,3) DO 30 III=1,LAST R=SQRT((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2) DL=R/EP1 IF(DL.GT.0.5*HX) DL=0.5*HX IF(DL.LT.DL0) DL=DL0 DL1=DL*AA(11) X2=X1+DL1*Q(M1)/B1 Y2=Y1+DL1*Q(M2)/B1 Z2=Z1+DL1*Q(M3)/B1 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 100 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 100 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 100 CALL QUANT1(X2,Y2,Z2,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 DO 32 II=1,1 X3=X1+0.5*DL1*(Q(M1)/B1+Q(M4)/B2) Y3=Y1+0.5*DL1*(Q(M2)/B1+Q(M5)/B2) Z3=Z1+0.5*DL1*(Q(M3)/B1+Q(M6)/B2) IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 100 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 100 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 100 CALL QUANT1(X3,Y3,Z3,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 32 CONTINUE C GX1=GX0+XDD1+AMP*(X3+AMPY*Y3*XC1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) CALL PLOT(GX1,GY1,2) X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE 30 CONTINUE 100 CONTINUE RETURN END SUBROUTINE AINTE2(IA,AA,F,P) DIMENSION IA(1),AA(1),F(1),P(1),Q(16) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO C AA 1 2 3 4 5 6 7 8 9 10 C TH0 ARU AR1 ARB XXL YYL ZZL B0 DL SIGN C 11 12 13 14 15 16 17 18 19 20 C DIRE BOUN GX0 GY0 GXL GYL GTH GXMI GXMA EP1 C P(1)=P(1) PI=3.1415926 NX=IA(1) NY=IA(2) MX=IA(3) MY=IA(4) MZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) TH0=AA(1)*PI/180.0 ARU=AA(2) AR1=AA(3) ARB=AA(4) XXL=AA(5) YYL=AA(6) ZZL=AA(7) B0=AA(8) GX0=AA(13) GY0=AA(14) GXL=AA(15) GYL=AA(16) GTH=AA(17)*PI/180.0 GXMI=AA(18) GXMA=AA(19) EP1=AA(20) C XX=999.0 XM=0.075 YM=0.075 X1=GX0 Y1=GY0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,MI,GXMI,GXMA) C XC1=COS(GTH) XS1=SIN(GTH) Y=GY0+0.5*GYL CALL PLOT(GX0,Y,3) X=GX0 Y=GY0+GYL CALL PLOT(X,Y,2) X=GX0+GXL CALL PLOT(X,Y,2) Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=GX0 CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1 Y=GY0+0.5*GYL*XS1+0.5*GYL CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1+GXL CALL PLOT(X,Y,2) X=GX0+GXL Y=GY0+0.5*GYL CALL PLOT(X,Y,2) C XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+GYL CALL PLOT(X,Y,3) Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=X+0.5*GYL*XC1 Y=Y+0.5*GYL*XS1 CALL PLOT(X,Y,2) C XL0=COS(TH0) XL01=XL0*1.1 XL=2.0*XL0 YL=XL NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 N1=NX2 N2=N1*NY2 AR2=AR1*AR1 DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) HX=XXL/FLOAT(MX+1) HY=YYL/FLOAT(MY+1) HZ=ZZL/FLOAT(MZ+1) XLMI=0.5*HX*FLOAT(2*NXP-MX-1) XLMA=0.5*HX*FLOAT(2*NXP+MX+1) DL0=HX*AA(9) DL=DL0 LAST=IFIX(XXL/DL) GYMI=-0.1*HY GYMA=0.5*(GXMA-GXMI) GZMI=-0.1*HZ GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 C JB1=NY2/2+1 DO 100 J=1,NY DO 100 I=1,NX1 X1=0.5*PI-FLOAT(J)*(0.5*PI-TH0)/FLOAT(NY) Y1=PI*FLOAT(I-1)/FLOAT(NX) R=COS(X1) X=R*COS(Y1) Y=R*SIN(Y1) IF(R.GT.XL01) GO TO 100 Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=ARB*ARB-X1*X1-Y1*Y1 IF(Z1.LE.0.0) GO TO 100 Z1=SQRT(Z1) CALL QUANT1(X1,Y1,Z1,HX,HY,HZ,1,MX,MY,MZ,NXP,AA,F,Q) B1=SQRT(Q(M1)*Q(M1)+Q(M2)*Q(M2)+Q(M3)*Q(M3)) IF(B1.LT.1.0E-20) GO TO 100 AA(11)=1.0 X2=Q(M1)*X1+Q(M2)*Y1+Q(M3)*Z1 IF(X2.LT.0.0) AA(11)=-1.0 C X2=0.0 Y2=0.0 Z2=0.0 X3=0.0 Y3=0.0 Z3=0.0 GX1=GX0+XDD1+AMP*(X1+Y1*XC1) GY1=GY0+0.5*GYL+AMP*(Z1+Y1*XS1) CALL PLOT(GX1,GY1,3) DO 30 III=1,LAST R=SQRT((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2) DL=R/EP1 IF(DL.GT.0.5*HX) DL=0.5*HX IF(DL.LT.DL0) DL=DL0 DL1=DL*AA(11) X2=X1+DL1*Q(M1)/B1 Y2=Y1+DL1*Q(M2)/B1 Z2=Z1+DL1*Q(M3)/B1 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 100 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 100 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 100 CALL QUANT1(X2,Y2,Z2,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 DO 32 II=1,1 X3=X1+0.5*DL1*(Q(M1)/B1+Q(M4)/B2) Y3=Y1+0.5*DL1*(Q(M2)/B1+Q(M5)/B2) Z3=Z1+0.5*DL1*(Q(M3)/B1+Q(M6)/B2) IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 100 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 100 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 100 CALL QUANT1(X3,Y3,Z3,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 32 CONTINUE C GX1=GX0+XDD1+AMP*(X3+Y3*XC1) GY1=GY0+0.5*GYL+AMP*(Z3+Y3*XS1) CALL PLOT(GX1,GY1,2) X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE 30 CONTINUE 100 CONTINUE RETURN END SUBROUTINE AINTE3(IA,AA,F,P) DIMENSION IA(1),AA(1),F(1),P(1),Q(16) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO C AA 1 2 3 4 5 6 7 8 9 10 C TH0 ARU AR1 ARB XXL YYL ZZL B0 DL SIGN C 11 12 13 14 15 16 17 18 19 20 C DIRE BOUN GX0 GY0 GXL GYL GTH GXMI GXMA EP1 C P(1)=P(1) PI=3.1415926 NX=IA(1) NY=IA(2) MX=IA(3) MY=IA(4) MZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) TH0=AA(1)*PI/180.0 ARU=AA(2) AR1=AA(3) ARB=AA(4) XXL=AA(5) YYL=AA(6) ZZL=AA(7) B0=AA(8) GX0=AA(13) GY0=AA(14) GXL=AA(15) GYL=AA(16) GTH=AA(17)*PI/180.0 GXMI=AA(18) GXMA=AA(19) EP1=AA(20) C XX=999.0 XM=0.075 YM=0.075 X1=GX0 Y1=GY0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,MI,GXMI,GXMA) C XC1=COS(GTH) XS1=SIN(GTH) X=GX0+0.5*GXL Y=GY0 CALL PLOT(X,Y,3) X=X-GXL*XC1 Y=Y-GXL*XS1 CALL PLOT(X,Y,2) X=X-0.5*GYL CALL PLOT(X,Y,2) X=GX0 Y=GY0 CALL PLOT(X,Y,2) X=GX0+0.5*GYL CALL PLOT(X,Y,2) Y=Y+0.5*GYL CALL PLOT(X,Y,2) X=X-GXL*XC1 Y=Y-GXL*XS1 CALL PLOT(X,Y,2) Y=Y-0.5*GYL CALL PLOT(X,Y,2) C XDD1=GXL*(GXMA-0.0)/(GXMA-GXMI) X=GX0-XDD1*XC1 Y=GY0-XDD1*XS1 CALL PLOT(X,Y,3) X=X+0.5*GYL CALL PLOT(X,Y,2) Y=Y+0.5*GYL CALL PLOT(X,Y,2) C XL0=COS(TH0) XL01=XL0*1.0 XL=2.0*XL0 YL=XL NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 N1=NX2 N2=N1*NY2 AR2=AR1*AR1 DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) HX=XXL/FLOAT(MX+1) HY=YYL/FLOAT(MY+1) HZ=ZZL/FLOAT(MZ+1) XLMI=0.5*HX*FLOAT(2*NXP-MX-1) XLMA=0.5*HX*FLOAT(2*NXP+MX+1) DL0=HX*AA(9) DL=DL0 LAST=IFIX(XXL/DL) GYMI=-0.1*HY GYMA=0.5*(GXMA-GXMI) GZMI=-0.1*HZ GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 C JB1=NY2/2+1 DO 100 J=JB1,NY2 DO 100 I=1,NX2 X=0.5*DX*(2*I-NX2-1) Y=0.5*DY*(2*J-NY2-1) R=SQRT(X*X+Y*Y) IF(R.GT.XL01) GO TO 100 Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=ARB*ARB-X1*X1-Y1*Y1 IF(Z1.LE.0.0) GO TO 100 Z1=SQRT(Z1) CALL QUANT1(X1,Y1,Z1,HX,HY,HZ,1,MX,MY,MZ,NXP,AA,F,Q) B1=SQRT(Q(M1)*Q(M1)+Q(M2)*Q(M2)+Q(M3)*Q(M3)) IF(B1.LT.1.0E-20) GO TO 100 AA(11)=1.0 X2=Q(M1)*X1+Q(M2)*Y1+Q(M3)*Z1 IF(X2.LT.0.0) AA(11)=-1.0 C X2=0.0 Y2=0.0 Z2=0.0 X3=0.0 Y3=0.0 Z3=0.0 GX1=GX0+0.5*GYL-XDD1*XC1+AMP*(-Y1+X1*XC1) GY1=GY0-XDD1*XS1+AMP*(Z1+X1*XS1) CALL PLOT(GX1,GY1,3) DO 30 III=1,LAST R=SQRT((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2) DL=R/EP1 IF(DL.GT.0.5*HX) DL=0.5*HX IF(DL.LT.DL0) DL=DL0 DL1=DL*AA(11) X2=X1+DL1*Q(M1)/B1 Y2=Y1+DL1*Q(M2)/B1 Z2=Z1+DL1*Q(M3)/B1 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 100 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 100 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 100 CALL QUANT1(X2,Y2,Z2,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 DO 32 II=1,1 X3=X1+0.5*DL1*(Q(M1)/B1+Q(M4)/B2) Y3=Y1+0.5*DL1*(Q(M2)/B1+Q(M5)/B2) Z3=Z1+0.5*DL1*(Q(M3)/B1+Q(M6)/B2) IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 100 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 100 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 100 CALL QUANT1(X3,Y3,Z3,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 32 CONTINUE C GX1=GX0+0.5*GYL-XDD1*XC1+AMP*(-Y3+X3*XC1) GY1=GY0-XDD1*XS1+AMP*(Z3+X3*XS1) CALL PLOT(GX1,GY1,2) X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE 30 CONTINUE 100 CONTINUE RETURN END SUBROUTINE AINTEG(IA,AA,F,P) DIMENSION IA(1),AA(1),F(1),P(1),Q(16) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO C AA 1 2 3 4 5 6 7 8 9 10 11 C TH0 ARU AR1 ARB XXL YYL ZZL B0 DL SIGN DIRE C PI=3.1415926 NX=IA(1) NY=IA(2) MX=IA(3) MY=IA(4) MZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) TH0=AA(1) ARU=AA(2) AR1=AA(3) ARB=AA(4) XXL=AA(5) YYL=AA(6) ZZL=AA(7) B0=AA(8) VXC=AA(22) C XL0=COS(PI*TH0/180.0) XL01=XL0*1.2 XL=2.0*XL0 YL=XL NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 N1=NX2 N2=N1*NY2 AR2=AR1*AR1 DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) HX=XXL/FLOAT(MX+1) HY=YYL/FLOAT(MY+1) HZ=ZZL/FLOAT(MZ+1) XLMI=0.5*HX*FLOAT(2*NXP-MX-1) XLMA=0.5*HX*FLOAT(2*NXP+MX+1) DL0=HX*AA(9) DL=DL0 EP1=1.0E-2 LAST=IFIX(XXL/DL) C JB1=NY2/2+1 DO 100 J=JB1,NY2 DO 100 I=1,NX2 X=0.5*DX*(2*I-NX2-1) Y=0.5*DY*(2*J-NY2-1) R=SQRT(X*X+Y*Y) IF(R.GT.XL01) GO TO 20 Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=SQRT(ARB*ARB-X1*X1-Y1*Y1) CALL QUANT1(X1,Y1,Z1,HX,HY,HZ,1,MX,MY,MZ,NXP,AA,F,Q) B1=SQRT(Q(6)*Q(6)+Q(7)*Q(7)+Q(8)*Q(8)) IF(B1.LT.1.0E-20) GO TO 20 AA(11)=1.0 X2=Q(6)*X1+Q(7)*Y1+Q(8)*Z1 IF(X2.LT.0.0) AA(11)=-1.0 C C FB1=0.0 FB5=0.0 BBL=0.0 X2=0.0 Y2=0.0 Z2=0.0 X3=0.0 Y3=0.0 Z3=0.0 DO 30 III=1,LAST R=SQRT((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2) DL=R/EP1 IF(DL.GT.HX) DL=HX IF(DL.LT.DL0) DL=DL0 DL1=DL*AA(11) X2=X1+DL1*Q(6)/B1 Y2=Y1+DL1*Q(7)/B1 Z2=Z1+DL1*Q(8)/B1 IF(X2.LT.XLMI.OR.X2.GT.XLMA) GO TO 20 IF(Y2.LT.0.0.OR.Y2.GT.YYL) GO TO 20 IF(Z2.LT.0.0.OR.Z2.GT.ZZL) GO TO 20 CALL QUANT1(X2,Y2,Z2,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(14)*Q(14)+Q(15)*Q(15)+Q(16)*Q(16)) IF(B2.LT.1.0E-20) GO TO 20 DO 32 II=1,1 X3=X1+0.5*DL1*(Q(6)/B1+Q(14)/B2) Y3=Y1+0.5*DL1*(Q(7)/B1+Q(15)/B2) Z3=Z1+0.5*DL1*(Q(8)/B1+Q(16)/B2) IF(X3.LT.XLMI.OR.X3.GT.XLMA) GO TO 20 IF(Y3.LT.0.0.OR.Y3.GT.YYL) GO TO 20 IF(Z3.LT.0.0.OR.Z3.GT.ZZL) GO TO 20 CALL QUANT1(X3,Y3,Z3,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(14)*Q(14)+Q(15)*Q(15)+Q(16)*Q(16)) IF(B2.LT.1.0E-20) GO TO 20 IF(Q(10).GT.VXC) GO TO 20 32 CONTINUE FB1=FB1+0.5*DL*(Q(1)/B1+Q(9)/B2) FB5=FB5+0.5*DL*(Q(5)/B1+Q(13)/B2) BBL=BBL+0.5*DL*(1.0/B1+1.0/B2) C BBL=BBL+0.5*DL X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE 30 CONTINUE 20 CONTINUE I1=J+NY2*(NX2-I)+N2*(MO-1) I2=NY2-J+1+NY2*(NX2-I)+N2*(MO-1) IF(BBL.LT.1.0E-20) P(I1)=0.0 IF(BBL.GE.1.0E-20) P(I1)=FB1/BBL P(I2)=P(I1)*AA(10) IF(BBL.LT.1.0E-20) P(I1+N2)=0.0 IF(BBL.GE.1.0E-20) P(I1+N2)=FB5/BBL P(I2+N2)=P(I1+N2)*AA(10) 100 CONTINUE RETURN END SUBROUTINE QUANT1(X,Y,Z,HX,HY,HZ,MM,NX,NY,NZ,NXP,AA,F,P) DIMENSION F(1),P(1),AA(1),B(3) C NX1=NX+1 NX2=NX+2 NY1=NY+1 NZ1=NZ+1 N1=NX+2 N2=N1*(NY+2) N3=N2*(NZ+2) N11=N1+1 N12=N2+1 N21=N2+N1 N22=N21+1 ARU=AA(2) AR1=AA(3) B0=AA(8) C XI=X/HX+0.5*FLOAT(NX2+1-2*NXP) XJ=Y/HY+1.5 XK=Z/HZ+1.5 I=IFIX(XI) J=IFIX(XJ) K=IFIX(XK) IF(I.LT.1) I=1 IF(J.LT.1) J=1 IF(K.LT.1) K=1 IF(I.GT.NX1) I=NX1 IF(J.GT.NY1) J=NY1 IF(K.GT.NZ1) K=NZ1 XI=XI-FLOAT(I) XJ=XJ-FLOAT(J) XK=XK-FLOAT(K) XI1=1.0-XI XJ1=1.0-XJ XK1=1.0-XK I1=I+N1*(J-1)+N2*(K-1)-N3 L1=8*(MM-1) DO 10 M=1,8 I1=I1+N3 I2=M+L1 P(I2)=XI*XJ*XK*F(I1+N22)+XI*XJ*XK1*F(I1+N11) & +XI*XJ1*XK*F(I1+N12)+XI1*XJ*XK*F(I1+N21) & +XI*XJ1*XK1*F(I1+1)+XI1*XJ*XK1*F(I1+N1) & +XI1*XJ1*XK*F(I1+N2)+XI1*XJ1*XK1*F(I1) 10 CONTINUE C AR2=AR1*AR1 RO2=X*X+Y*Y+Z*Z RO5=RO2*RO2*SQRT(RO2) X2=0.0 IF(RO2.GT.AR2) X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) B(1)=-3.0*B0*X*Z/RO5 B(2)=-3.0*B0*Y*Z/RO5 B(3)=B0*(RO2-3.0*Z*Z)/RO5 L1=5+8*(MM-1) DO 20 M=1,3 I1=M+L1 P(I1)=P(I1)*X2+B(M)*(1.0-X2) 20 CONTINUE RETURN END SUBROUTINE ALONGB(IA,AA,F,P) DIMENSION IA(1),AA(1),F(1),P(1) P(1)=P(1) AA(1)=AA(1) NX=IA(3) NY=IA(4) NZ=IA(5) MI=IA(7) MO=IA(8) NX2=NX+2 NY2=NY+2 NZ2=NZ+2 N1=NX+2 N2=N1*(NY+2) N3=N2*(NZ+2) N32=N3*2 DO 10 K=1,NZ2 DO 10 J=1,NY2 L1=N1*(J-1)+N2*(K-1)+N3*(MI-1) L2=N1*(J-1)+N2*(K-1)+N3*(MO-1) L3=N1*(J-1)+N2*(K-1)+N3*5 DO 10 I=1,NX2 I6=I+L3 I7=I6+N3 I8=I7+N3 II=I+L1 IO=I+L2 X=SQRT(F(I6)*F(I6)+F(I7)*F(I7)+F(I8)*F(I8)) IF(X.LT.1.0E-20) F(IO)=0.0 IF(X.GE.1.0E-20) F(IO)=(F(I6)*F(II)+F(I7)*F(II+N3) & +F(I8)*F(II+N32))/X 10 CONTINUE RETURN END SUBROUTINE PARAB1(IA,AA,F,P) DIMENSION IA(1),AA(1),F(1),P(1) P(1)=P(1) AA(1)=AA(1) NX=IA(3) NY=IA(4) NZ=IA(5) MI=IA(7) MO=IA(8) NX2=NX+2 NY2=NY+2 NZ2=NZ+2 N1=NX+2 N2=N1*(NY+2) N3=N2*(NZ+2) N32=N3*2 DO 10 K=1,NZ2 DO 10 J=1,NY2 L1=N1*(J-1)+N2*(K-1)+N3*(MI-1) L2=N1*(J-1)+N2*(K-1)+N3*(MO-1) L3=N1*(J-1)+N2*(K-1)+N3*5 DO 10 I=1,NX2 I6=I+L3 I7=I6+N3 I8=I7+N3 II=I+L1 IO=I+L2 X=F(I6)*F(I6)+F(I7)*F(I7)+F(I8)*F(I8) IF(X.LT.1.0E-20) X1=0.0 IF(X.GE.1.0E-20) X1=(F(I6)*F(II)+F(I7)*F(II+N3) & +F(I8)*F(II+N32))/X F(II)=X1*F(I6) F(II+N3)=X1*F(I7) F(II+N32)=X1*F(I8) 10 CONTINUE RETURN END SUBROUTINE ROTAT(IA,AA,F,P) DIMENSION IA(1),AA(1),F(1),P(1) C NX=IA(3) NY=IA(4) NZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) ARU=AA(2) AR1=AA(3) XXL=AA(5) YYL=AA(6) ZZL=AA(7) NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 NZ1=NZ+1 NZ2=NZ+2 N1=NX2 N2=N1*NY2 N3=N2*NZ2 N11=N1+1 N12=N2+1 N21=N1+N2 N22=N21+1 AR2=AR1*AR1 HX=XXL/FLOAT(NX1) HY=YYL/FLOAT(NY1) HZ=ZZL/FLOAT(NZ1) C DO 10 K=2,NZ1 DO 12 M=1,3 DO 12 J=1,NY1 L1=N1*(J-1)+N2*(K-1)+N3*(M+MI-2) L2=N1*(J-1)+N2+2*N2*(M-1) DO 12 I=1,NX1 I1=I+L1 I2=I1-N2 J1=I+L2 IF(K.EQ.2) P(J1-N2)=0.125*(F(I2)+F(I2+1)+F(I2+N1)+F(I2+N11) & +F(I2+N2)+F(I2+N12)+F(I2+N21)+F(I2+N22)) IF(K.GE.3) P(J1-N2)=P(J1) P(J1)=0.125*(F(I1)+F(I1+1)+F(I1+N1)+F(I1+N11) & +F(I1+N2)+F(I1+N12)+F(I1+N21)+F(I1+N22)) 12 CONTINUE C C CURRENT DO 10 J=2,NY1 L1=N1*(J-1)+N2*(K-1)+N3*(MO-1) L2=N1*(J-1)+N2 DO 10 I=2,NX1 I6=I+L2 I7=I6+2*N2 I8=I7+2*N2 J2=I+L1 J3=J2+N3 J4=J3+N3 F(J2)=0.25*((P(I8)-P(I8-N1)+P(I8-N2)-P(I8-N21) & +P(I8-1)-P(I8-N11)+P(I8-N12)-P(I8-N22))/HY & -(P(I7)+P(I7-N1)-P(I7-N2)-P(I7-N21) & +P(I7-1)+P(I7-N11)-P(I7-N12)-P(I7-N22))/HZ) F(J3)=0.25*((P(I6)+P(I6-N1)-P(I6-N2)-P(I6-N21) & +P(I6-1)+P(I6-N11)-P(I6-N12)-P(I6-N22))/HZ & -(P(I8)+P(I8-N1)+P(I8-N2)+P(I8-N21) & -P(I8-1)-P(I8-N11)-P(I8-N12)-P(I8-N22))/HX) F(J4)=0.25*((P(I7)+P(I7-N1)+P(I7-N2)+P(I7-N21) & -P(I7-1)-P(I7-N11)-P(I7-N12)-P(I7-N22))/HX & -(P(I6)-P(I6-N1)+P(I6-N2)-P(I6-N21) & +P(I6-1)-P(I6-N11)+P(I6-N12)-P(I6-N22))/HY) 10 CONTINUE C C BOUNDARY CONDITION AT NZ=1 AND NZ=NZ2 C BOUNDARY CONDITION AT NY=1 AND NY=NY2 IF(MI.EQ.2) GO TO 20 DO 30 M=1,3 DO 31 J=2,NY1 DO 31 I=2,NX1 I1=I+N1*(J-1)+N3*(M+MO-2) I2=I1+N2*NZ1 I3=N2*2 IF(I.EQ.2) F(I2)=AA(12)*F(I2-N2-1) IF(I.GT.2) F(I2)=AA(12)*(4.0*F(I2-N2-1)-F(I2-I3-2))/3.0 F(I1)=F(I1+N2) C IF(M.EQ.1) F(I1)=-F(I1) IF(M.EQ.3) F(I1)=-F(I1) 31 CONTINUE DO 32 K=1,NZ2 DO 32 I=2,NX1 I1=I+N2*(K-1)+N3*(M+MO-2) I2=I1+N1*NY1 I3=N1*2 IF(I.EQ.2) F(I2)=AA(12)*F(I2-N1-1) IF(I.GT.2) F(I2)=AA(12)*(4.0*F(I2-N1-1)-F(I2-I3-2))/3.0 F(I1)=F(I1+N1) IF(M.EQ.1) F(I1)=-F(I1) IF(M.EQ.3) F(I1)=-F(I1) 32 CONTINUE C BOUNDARY CONDITION AT NX=NX2 DO 33 K=1,NZ2 DO 33 J=1,NY2 I1=1+N1*(J-1)+N2*(K-1)+N3*(M+MO-2) I2=I1+NX1 F(I1)=0.0 F(I2)=(4.0*F(I2-1)-F(I2-2))/3.0 33 CONTINUE 30 CONTINUE GO TO 60 C 20 CONTINUE DO 50 M=1,3 DO 51 J=2,NY1 DO 51 I=2,NX1 I1=I+N1*(J-1)+N3*(M+MO-2) I2=I1+N2*NZ1 I3=N2*2 IF(I.EQ.2) F(I2)=AA(12)*F(I2-N2-1) IF(I.GT.2) F(I2)=AA(12)*(4.0*F(I2-N2-1)-F(I2-I3-2))/3.0 F(I1)=F(I1+N2) IF(M.EQ.1) F(I1)=-F(I1) IF(M.EQ.2) F(I1)=-F(I1) 51 CONTINUE DO 52 K=1,NZ2 DO 52 I=2,NX1 I1=I+N2*(K-1)+N3*(M+MO-2) I2=I1+N1*NY1 I3=N1*2 IF(I.EQ.2) F(I2)=AA(12)*F(I2-N1-1) IF(I.GT.2) F(I2)=AA(12)*(4.0*F(I2-N1-1)-F(I2-I3-2))/3.0 F(I1)=F(I1+N1) IF(M.EQ.1) F(I1)=-F(I1) IF(M.EQ.3) F(I1)=-F(I1) 52 CONTINUE C BOUNDARY CONDITION AT NX=NX2 DO 53 K=1,NZ2 DO 53 J=1,NY2 I1=1+N1*(J-1)+N2*(K-1)+N3*(M+MO-2) I2=I1+NX1 F(I1)=0.0 F(I2)=(4.0*F(I2-1)-F(I2-2))/3.0 53 CONTINUE 50 CONTINUE C 60 CONTINUE DO 64 K=1,NZ2 DO 64 J=1,NY2 DO 64 I=1,NX2 X=0.5*HX*FLOAT(2*I-NX1-2+2*NXP) Y=0.5*HY*FLOAT(2*J-3) Z=0.5*HZ*FLOAT(2*K-3) RO2=X*X+Y*Y+Z*Z I2=I+N1*(J-1)+N2*(K-1)+N3*(MO-1) I3=I2+N3 I4=I3+N3 X2=RO2/AR2-1.0 IF(X2.LT.0.0) X2=0.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) F(I2)=F(I2)*X2 F(I3)=F(I3)*X2 F(I4)=F(I4)*X2 64 CONTINUE RETURN END SUBROUTINE PHAS2(III,NX,NZ,F,U) DIMENSION F(1),U(1) MX=2*NX+3 MY=2*NZ+3 NX1=NX+1 NX2=NX+2 NZ1=NZ+1 NZ2=NZ+2 N1=NX2*NZ2 N2=MX*MY MX3=MX*3 DO 10 I=1,N2 U(I)=0.0 10 CONTINUE DO 20 J=1,NZ2 DO 20 I=1,NX2 I1=2*I-1+MX*(2*J-2) I2=I+NX2*(J-1)+N1*(III-1) U(I1)=F(I2) 20 CONTINUE DO 30 J=1,NZ2 J1=2*J-1 DO 30 I=1,NX1 I1=2*I K=I1+MX*(J1-1) U(K)=0.5*(U(K-1)+U(K+1)) 30 CONTINUE DO 40 I=1,MX DO 40 J=1,NZ1 J1=2*J K=I+MX*(J1-1) U(K)=0.5*(U(K-MX)+U(K+MX)) 40 CONTINUE RETURN END SUBROUTINE DIFF2(M,MX,MY,U,V) DIMENSION U(1),V(1) C IF(M.EQ.99) GO TO 99 ALP=1.0/8.0 MX1=MX-1 MY1=MY-1 N1=MX*MY DO 100 II=1,M DO 12 I=1,N1 V(I)=U(I) 12 CONTINUE DO 14 J=2,MY1 DO 14 I=2,MX1 I1=I+MX*(J-1) U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+V(I1-MX)+V(I1+MX)) 14 CONTINUE MY2=MX*(MY-2) C DO 16 I=2,MX1 I1=I U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1+MX)) I1=I+MX*(MY-1) U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1-MX)) 16 CONTINUE DO 18 J=2,MY1 I1=1+MX*(J-1) U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1+1)+V(I1-MX)+V(I1+MX)) I1=MX+MX*(J-1) U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1-1)+V(I1-MX)+V(I1+MX)) 18 CONTINUE C DO 20 I=1,MX I1=I+MX*MY1/2 U(I1)=V(I1) 20 CONTINUE C 100 CONTINUE 99 CONTINUE RETURN END SUBROUTINE PHASS(III,NR,NZ,F,U) DIMENSION F(1),U(1) NX=2*NR+3 NY=2*NZ+3 NR1=NR+1 NR2=NR+2 NZ1=NZ+1 NZ2=NZ+2 N1=NR2*NZ2 N2=NX*NY NX3=NX*3 DO 10 I=1,N2 U(I)=0.0 10 CONTINUE DO 20 J=1,NZ2 DO 20 I=1,NR2 I1=2*I-1+NX*(2*J-2) I2=I+NR2*(J-1)+N1*(III-1) U(I1)=F(I2) 20 CONTINUE DO 30 J=1,NZ2 J1=2*J-1 I1=2 K=I1+NX*(J1-1) U(K)=(3.0*U(K-1)+6.0*U(K+1)-U(K+3))/8.0 I1=2*NR1 K=I1+NX*(J1-1) U(K)=(-U(K-3)+6.0*U(K-1)+3.0*U(K+1))/8.0 DO 30 I=2,NR I1=2*I K=I1+NX*(J1-1) U(K)=(-U(K-3)+9.0*U(K-1)+9.0*U(K+1)-U(K+3))/16.0 30 CONTINUE DO 40 I=1,NX J1=2 K=I+NX*(J1-1) U(K)=(3.0*U(K-NX)+6.0*U(K+NX)-U(K+NX3))/8.0 J1=2*NZ1 K=I+NX*(J1-1) U(K)=(-U(K-NX3)+6.0*U(K-NX)+3.0*U(K+NX))/8.0 DO 40 J=2,NZ J1=2*J K=I+NX*(J1-1) U(K)=(-U(K-NX3)+9.0*U(K-NX)+9.0*U(K+NX)-U(K+NX3))/16.0 40 CONTINUE RETURN END SUBROUTINE DIFFU(M,NX,NY,U,V) DIMENSION U(1),V(1) C IF(M.EQ.99) GO TO 99 ALP=1.0/8.0 NX1=NX-1 NY1=NY-1 N1=NX*NY DO 100 II=1,M DO 12 I=1,N1 V(I)=U(I) 12 CONTINUE DO 14 J=2,NY1 DO 14 I=2,NX1 I1=I+NX*(J-1) U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+V(I1-NX)+V(I1+NX)) 14 CONTINUE NY2=NX*(NY-2) C DO 16 I=2,NX1 I1=I U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1+NX)) I1=I+NX*(NY-1) U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1-NX)) 16 CONTINUE DO 18 J=2,NY1 I1=1+NX*(J-1) C U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1+1)+V(I1-NX)+V(I1+NX)) I1=NX+NX*(J-1) U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1-1)+V(I1-NX)+V(I1+NX)) 18 CONTINUE C 100 CONTINUE 99 CONTINUE RETURN END SUBROUTINE GAROHD(X0,Y0,X1,Y1,B,CR) C=CR*B EP0=1.0E-10 EP1=1.0E-3 EP2=1.0E3 IF(ABS(X1-X0).LT.EP0) GO TO 2 IF(ABS(Y1-Y0).LT.EP0) GO TO 3 AA=(Y1-Y0)/(X1-X0) IF(ABS(AA).GT.EP2) GO TO 2 IF(ABS(AA).LT.EP1) GO TO 3 AA2=AA*AA A1=B/SQRT(1.0+AA2) A2=C/SQRT(1.0+1.0/AA2) C XH=X1-A1 IF(X1.LT.X0) XH=X1+A1 X2=XH+A2 X3=XH-A2 YH=Y1+AA*(XH-X1) Y2=YH-(X2-XH)/AA Y3=YH-(X3-XH)/AA GO TO 1 2 CONTINUE X2=X1-C X3=X1+C Y2=Y1-B IF(Y1.LT.Y0) Y2=Y1+B Y3=Y2 GO TO 1 3 CONTINUE Y2=Y1-C Y3=Y1+C X2=X1-B IF(X1.LT.X0) X2=X1+B X3=X2 1 CONTINUE C CALL PLOT(X0,Y0,3) CALL PLOT(X1,Y1,2) CALL PLOT(X2,Y2,3) CALL PLOT(X1,Y1,2) CALL PLOT(X3,Y3,2) RETURN END SUBROUTINE GCIRCL(X,Y,R) DIMENSION SC(10) DATA SC/0.0,0.174,0.342,0.5,0.643,0.766,0.866,0.94,0.985,1.0/ X1=X+R CALL PLOT(X1,Y,3) DO 10 I=1,9 X1=X+SC(10-I)*R Y1=Y+SC(I+1)*R CALL PLOT(X1,Y1,2) 10 CONTINUE DO 20 I=10,18 X1=X-SC(I-8)*R Y1=Y+SC(19-I)*R CALL PLOT(X1,Y1,2) 20 CONTINUE DO 30 I=19,27 X1=X-SC(28-I)*R Y1=Y-SC(I-17)*R CALL PLOT(X1,Y1,2) 30 CONTINUE DO 40 I=28,36 X1=X+SC(I-26)*R Y1=Y-SC(37-I)*R CALL PLOT(X1,Y1,2) 40 CONTINUE RETURN END SUBROUTINE GCIRC1(X,Y,R,N) DIMENSION SC(90) PI=3.1415926 N1=N-1 DO 5 I=1,N X1=0.5*PI*FLOAT(I-1)/FLOAT(N1) 5 SC(I)=SIN(X1) X1=X+R CALL PLOT(X1,Y,3) DO 10 I=1,N1 X1=X+SC(N-I)*R Y1=Y+SC(I+1)*R CALL PLOT(X1,Y1,2) 10 CONTINUE DO 20 I=1,N1 X1=X-SC(I+1)*R Y1=Y+SC(N-I)*R CALL PLOT(X1,Y1,2) 20 CONTINUE DO 30 I=1,N1 X1=X-SC(N-I)*R Y1=Y-SC(I+1)*R CALL PLOT(X1,Y1,2) 30 CONTINUE DO 40 I=1,N1 X1=X+SC(I+1)*R Y1=Y-SC(N-I)*R CALL PLOT(X1,Y1,2) 40 CONTINUE RETURN END SUBROUTINE GRECT(X,Y,A,CR) ACR=A*CR X1=X-A Y1=Y-ACR X2=X+A Y2=Y+ACR CALL PLOT(X1,Y1,3) CALL PLOT(X1,Y2,2) CALL PLOT(X2,Y2,2) CALL PLOT(X2,Y1,2) CALL PLOT(X1,Y1,2) RETURN END SUBROUTINE GREC1(X,Y,A,CR,L) DA=A/FLOAT(L) DO 10 I=1,L A1=DA*FLOAT(I) ACR=A1*CR X1=X-A1 Y1=Y-ACR X2=X+A1 Y2=Y+ACR CALL PLOT(X1,Y1,3) CALL PLOT(X1,Y2,2) CALL PLOT(X2,Y2,2) CALL PLOT(X2,Y1,2) CALL PLOT(X1,Y1,2) 10 CONTINUE RETURN END C SCALL1 SUBROUTINE SCALL1(JXG,JYG,NX2,YMIN,YMAX,P) DIMENSION P(1) YMIN=P(1) YMAX=P(1) DO 532 J=1,JYG DO 532 I=1,JXG I1=I+NX2*(J-1) Y=P(I1) YMIN=AMIN1(YMIN,Y) YMAX=AMAX1(YMAX,Y) 532 CONTINUE RETURN END C SUBROUTINE CONT11(NX,NY,U) DIMENSION U(1) DO 40 J=1,NY DO 40 I=1,NX I1=I+NX*(J-1) U(I1+NX*NY)=U(I1) 40 CONTINUE DO 42 J=1,NY DO 42 I=1,NX I1=I+1+(NX+2)*J I2=I+NX*(J-1)+NX*NY U(I1)=U(I2) 42 CONTINUE NY1=NY+1 DO 44 J=2,NY1 I1=1+(NX+2)*(J-1) I2=I1+NX+1 U(I1)=U(I1+NX) U(I2)=U(I2-NX) 44 CONTINUE NX2=NX+2 DO 46 I=1,NX2 I1=I I2=I+(NX+2)*(NY+1) I3=(NX+2) U(I1)=U(I1+I3) U(I2)=U(I2-I3) 46 CONTINUE RETURN END C /GR/PHASP SUBROUTINE PHASP(III,NX,NY,NXX,VE,F,U) DIMENSION F(1),U(1) PI=3.1415926 PI1=1.0/SQRT(2.0*PI) N1=NXX+2 NY2=NY/2 INX=NXX/NX N2=NX*NY DO 10 I=1,N2 U(I)=0.0 10 CONTINUE C DO 20 I=1,NX DO 20 J=1,NY V=FLOAT(J-NY2)*VE/FLOAT(NY2) DO 20 II=1,2 I1=I*INX+3*N1*(II-1)+6*N1*(III-1) XN=F(I1) XV=F(I1+N1) XVT2=F(I1+2*N1) XVT=SQRT(XVT2) X=0.5*(V-XV)*(V-XV)/XVT2 A=0.0 IF(X.LT.5.0) A=PI1*XN*EXP(-X)/XVT I2=I+NX*(J-1) U(I2)=U(I2)+A 20 CONTINUE RETURN END SUBROUTINE GRAP6G(IA,LB,LL,NX,NXX,X0,Y0,XL,YL,U,VMIN,VMAX,KP) C DIMENSION U(1) XX=999.0 XM=0.075 YM=0.075 C NX1=NX+1 NX2=NX+2 N1=NXX+2 C INX=NXX/NX L1=LB+LL-1 DO 8 J=LB,L1 DO 8 I=1,NX I1=I+N1*(J-1) I2=I*INX+N1*(J-1) U(I1)=U(I2) 8 CONTINUE DX=XL/FLOAT(NX1) I1=1+N1*(LB-1) UMAX=U(I1)+1.0E-6 UMIN=U(I1)-1.0E-6 DO 10 J=LB,L1 DO 10 I=1,NX I1=I+N1*(J-1) UMAX=AMAX1(UMAX,U(I1)) UMIN=AMIN1(UMIN,U(I1)) 10 CONTINUE IF(KP.EQ.0) GO TO 16 IF(KP.EQ.1) GO TO 12 IF(KP.EQ.-1) GO TO 14 IF(KP.EQ.5) GO TO 15 UMIN=VMIN UMAX=VMAX GO TO 16 12 UMAX=VMAX GO TO 16 14 UMIN=VMIN GO TO 16 15 UMIN=-UMIN UMAX=AMAX1(UMAX,UMIN) UMIN=-UMAX 16 CONTINUE C DO 18 J=LB,L1 DO 18 I=1,NX I1=I+N1*(J-1) IF(U(I1).GT.UMAX) U(I1)=UMAX IF(U(I1).LT.UMIN) U(I1)=UMIN 18 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,1,IA,UMIN,UMAX) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) C DO 20 JJ=LB,L1 X=X0 I1=NX+N1*(JJ-1) Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,3) DO 30 I=1,NX I1=I+N1*(JJ-1) X=X+DX Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,2) 30 CONTINUE I1=1+N1*(JJ-1) X=X+DX Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,2) 20 CONTINUE RETURN END SUBROUTINE GRAP4O(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX, & KP,LANK,V0,TH0,DTH,IPEN,AFAC) C DIMENSION U(1) TH0=TH0 DTH=DTH PI=3.1415926 XX=999.0 XM=0.075 YM=0.075 NXY=NX*NY NX1=NX+1 NY1=NY+1 NX2=NX+2 NY2=NY+2 N1=NX2 N2=N1*NY2 IHI=3 ILO=2 CALL NEWPE1(IPEN) C DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 10 J=1,NY2 DO 10 I=1,NX2 Y=0.5*HY*FLOAT(2*J-NY-3) X=0.5*HX*FLOAT(2*I-NX-3+2*NXP) R=X*X+Y*Y IF(R.LT.AR2) GO TO 10 I1=I+NX2*(J-1) X1=U(I1) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 10 CONTINUE IF(KP.EQ.0) GO TO 16 IF(KP.EQ.1) GO TO 12 IF(KP.EQ.-1) GO TO 14 IF(KP.EQ.5) GO TO 15 UMIN=VMIN UMAX=VMAX GO TO 16 12 UMAX=VMAX GO TO 16 14 UMIN=VMIN GO TO 16 15 CONTINUE X=ABS(UMAX-V0) Y=ABS(V0-UMIN) X=AMAX1(X,Y) UMAX=V0+X UMIN=V0-X 16 CONTINUE C DO 18 I=1,N2 IF(U(I).GT.UMAX) U(I)=UMAX IF(U(I).LT.UMIN) U(I)=UMIN 18 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,KP,UMIN,UMAX) C CALL NEWPE1(ILO) CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) Y=Y0+0.5*YL CALL PLOT(X,Y,3) X=X0+XL CALL PLOT(X,Y,2) X=X0+XL*(0.5-FLOAT(NXP)/FLOAT(NX1)) Y=Y0 CALL PLOT(X,Y,3) Y=Y+YL CALL PLOT(X,Y,2) C ULD=FLOAT(LANK+1)/(UMAX-UMIN) U00=(V0-UMIN)*ULD EP1=1.0E-8 XO0=X0+0.5*XL YO0=Y0+0.5*YL RO0=(0.55*XL)**2 C DO 100 J=1,NY1 DO 100 I=1,NX1 C I1=I+N1*(J-1) I2=I1+1 I3=I1+1+N1 I4=I1+N1 C U1=(U(I1)-UMIN)*ULD U2=(U(I2)-UMIN)*ULD U3=(U(I3)-UMIN)*ULD U4=(U(I4)-UMIN)*ULD X1=X0+DX*FLOAT(I-1) X2=X1+DX Y1=Y0+DY*FLOAT(J-1) Y2=Y1+DY U0=0.25*(U1+U2+U3+U4) X=U00+AFAC*(UMAX-UMIN)*ULD IF(U0.LT.X) GO TO 50 CALL NEWPE1(1) CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) CALL PLOT(X1,Y2,3) CALL PLOT(X2,Y1,2) 50 CONTINUE C K1=IFIX(U1) K2=IFIX(U2) K3=IFIX(U3) K4=IFIX(U4) C J1=IABS(K2-K1) J2=IABS(K3-K2) J3=IABS(K4-K3) J4=IABS(K1-K4) C IF(J1.EQ.0) GO TO 28 DO 20 I1=1,J1 U0=FLOAT(I1)+FLOAT(MIN0(K1,K2)) IF(U0.GE.U00) CALL NEWPE1(IHI) IF(U0.LT.U00) CALL NEWPE1(ILO) IF(ABS(U2-U1).LT.EP1) GO TO 20 X1=DX*((U0-U1)/(U2-U1)+FLOAT(I-1)) Y1=DY*FLOAT(J-1) X1=X0+X1 Y1=Y0+Y1 C IF((U3-U0)*(U2-U0).GT.0.0) GO TO 22 IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 22 IF(ABS(U3-U2).LT.EP1) GO TO 22 X2=DX*FLOAT(I) Y2=DY*((U0-U2)/(U3-U2)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 22 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 22 CONTINUE C IF((U4-U0)*(U3-U0).GT.0.0) GO TO 24 IF((U1-U0)*(U3-U0).GT.0.0) GO TO 24 IF((U2-U0)*(U4-U0).GT.0.0) GO TO 24 IF(ABS(U4-U3).LT.EP1) GO TO 24 X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y2=DY*FLOAT(J) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 24 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 24 CONTINUE C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 26 IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 26 IF(ABS(U1-U4).LT.EP1) GO TO 26 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 26 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 26 CONTINUE 20 CONTINUE 28 CONTINUE C IF(J2.EQ.0) GO TO 38 DO 30 I1=1,J2 IF(U0.GE.U00) CALL NEWPE1(IHI) IF(U0.LT.U00) CALL NEWPE1(ILO) U0=FLOAT(I1)+FLOAT(MIN0(K2,K3)) IF(ABS(U3-U2).LT.EP1) GO TO 30 X1=DX*FLOAT(I) Y1=DY*((U0-U2)/(U3-U2)+FLOAT(J-1)) X1=X0+X1 Y1=Y0+Y1 C IF((U4-U0)*(U3-U0).GT.0.0) GO TO 32 IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 32 IF(ABS(U4-U3).LT.EP1) GO TO 32 X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y2=DY*FLOAT(J) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 32 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 32 CONTINUE C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 34 IF((U1-U0)*(U3-U0).GT.0.0) GO TO 34 IF((U2-U0)*(U4-U0).GT.0.0) GO TO 34 IF(ABS(U1-U4).LT.EP1) GO TO 34 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 34 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 34 CONTINUE 30 CONTINUE 38 CONTINUE C IF(J3.EQ.0) GO TO 48 DO 40 I1=1,J3 U0=FLOAT(I1)+FLOAT(MIN0(K3,K4)) IF(U0.GE.U00) CALL NEWPE1(IHI) IF(U0.LT.U00) CALL NEWPE1(ILO) IF(ABS(U4-U3).LT.EP1) GO TO 40 X1=DX*((U0-U4)/(U3-U4)+FLOAT(I-1)) Y1=DY*FLOAT(J) X1=X0+X1 Y1=Y0+Y1 C IF((U1-U0)*(U4-U0).GT.0.0) GO TO 42 IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 42 IF(ABS(U1-U4).LT.EP1) GO TO 42 X2=DX*FLOAT(I-1) Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1)) X2=X0+X2 Y2=Y0+Y2 R1=(X1-XO0)**2+(Y1-YO0)**2 R2=(X2-XO0)**2+(Y2-YO0)**2 IF(R1.GT.RO0.AND.R2.GT.RO0) GO TO 42 CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) 42 CONTINUE 40 CONTINUE 48 CONTINUE C 100 CONTINUE CALL NEWPE1(IPEN) RETURN END SUBROUTINE GRAP1N(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN, & VMAX,URMI0,AD,TH1,KP,KPS,IAR,XEP) DIMENSION U(1) C XX=999.0 XM=0.075 YM=0.075 NXY=NX*NY NL=NXY*2 DX=XL/FLOAT(NX+1) DY=YL/FLOAT(NY+1) AR=DX*0.40 CR=0.40 UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 8 J=1,NY DO 8 I=1,NX Y=0.5*HY*FLOAT(2*J-NY-1) X=0.5*HX*FLOAT(2*I-NX-1+2*NXP) R=X*X+Y*Y IF(R.LT.AR2) GO TO 8 I1=I+NX*(J-1) I2=I1+NXY X1=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 8 CONTINUE IF(KP.EQ.0) GO TO 126 IF(KP.EQ.1) GO TO 122 IF(KP.EQ.-1) GO TO 124 UMIN=VMIN UMAX=VMAX GO TO 126 122 UMAX=VMAX GO TO 126 124 UMIN=VMIN 126 CONTINUE C UMIN=-UMIN UMAX=AMAX1(UMAX,UMIN) UMIN=-UMAX URMI=URMI0*UMAX C DO 18 I=1,NXY I2=I+NXY R=SQRT(U(I)*U(I)+U(I2)*U(I2)) IF(R.GT.UMAX) U(I)=U(I)*UMAX/R IF(R.GT.UMAX) U(I2)=U(I2)*UMAX/R 18 CONTINUE C IF(KPS.EQ.0) GO TO 20 X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,KP,UMIN,UMAX) C XC1=COS(TH1) XS1=SIN(TH1) Y=Y0+0.5*YL CALL PLOT(X0,Y,3) X=X0 Y=Y0+YL CALL PLOT(X,Y,2) X=X0+XL CALL PLOT(X,Y,2) Y=Y0+0.5*YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) X=X0+0.5*YL*XC1 Y=Y0+0.5*YL*XS1+0.5*YL CALL PLOT(X,Y,2) X=X0+0.5*YL*XC1+XL CALL PLOT(X,Y,2) X=X0+XL Y=Y0+0.5*YL CALL PLOT(X,Y,2) X=X0+XEP*XL Y=Y0+YL CALL PLOT(X,Y,3) Y=Y0+0.5*YL CALL PLOT(X,Y,2) X=X+0.5*YL*XC1 Y=Y0+0.5*YL*XS1+0.5*YL CALL PLOT(X,Y,2) 20 CONTINUE C NY1=NY/2 DO 10 J=1,NY1 Y1=Y0+DY*FLOAT(2*J-1) DO 12 I=1,NX X1=X0+DX*FLOAT(I) I1=I+NX*(2*J-2) I2=I1+NXY R=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) IF(R.LT.URMI) GO TO 12 UX=AD*DX*U(I1)/UMAX UY=AD*DY*U(I2)/UMAX X2=X1+UX Y2=Y1+UY XX1=SQRT((UX*UX+UY*UY)/(DX*DX+DY*DY)) AR1=AR IF(IAR.EQ.1) AR1=AR*XX1*1.5 CALL GAROHD(X1,Y1,X2,Y2,AR1,CR) 12 CONTINUE Y1=Y1+DY DO 14 I=1,NX X1=X0+DX*FLOAT(NX+1-I) I1=NX+1-I+NX*(2*J-1) I2=I1+NXY R=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) IF(R.LT.URMI) GO TO 14 UX=AD*DX*U(I1)/UMAX UY=AD*DY*U(I2)/UMAX X2=X1+UX Y2=Y1+UY XX1=SQRT((UX*UX+UY*UY)/(DX*DX+DY*DY)) AR1=AR IF(IAR.EQ.1) AR1=AR*XX1*1.5 CALL GAROHD(X1,Y1,X2,Y2,AR1,CR) 14 CONTINUE 10 CONTINUE RETURN END SUBROUTINE GRAP1M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN, & VMAX,URMI0,KP,IAR,XEP) DIMENSION U(1) C XX=999.0 XM=0.075 YM=0.075 NXY=NX*NY NL=NXY*2 AD=1.0 DX=XL/FLOAT(NX+1) DY=YL/FLOAT(NY+1) AR=DX*0.40 CR=0.40 UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 8 J=1,NY DO 8 I=1,NX Y=0.5*HY*FLOAT(2*J-NY-1) X=0.5*HX*FLOAT(2*I-NX-1+2*NXP) R=X*X+Y*Y IF(R.LT.AR2) GO TO 8 I1=I+NX*(J-1) I2=I1+NXY X1=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 8 CONTINUE IF(KP.EQ.0) GO TO 126 IF(KP.EQ.1) GO TO 122 IF(KP.EQ.-1) GO TO 124 UMIN=VMIN UMAX=VMAX GO TO 126 122 UMAX=VMAX GO TO 126 124 UMIN=VMIN 126 CONTINUE C UMIN=-UMIN UMAX=AMAX1(UMAX,UMIN) UMIN=-UMAX URMI=URMI0*UMAX C DO 18 I=1,NXY I2=I+NXY R=SQRT(U(I)*U(I)+U(I2)*U(I2)) IF(R.GT.UMAX) U(I)=U(I)*UMAX/R IF(R.GT.UMAX) U(I2)=U(I2)*UMAX/R 18 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,KP,UMIN,UMAX) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) Y=Y0+0.5*YL CALL PLOT(X,Y,3) X=X0+XL CALL PLOT(X,Y,2) X=X0+XEP*XL Y=Y0 CALL PLOT(X,Y,3) Y=Y0+YL CALL PLOT(X,Y,2) C NY1=NY/2 DO 10 J=1,NY1 Y1=Y0+DY*FLOAT(2*J-1) DO 12 I=1,NX X1=X0+DX*FLOAT(I) I1=I+NX*(2*J-2) I2=I1+NXY R=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) IF(R.LT.URMI) GO TO 12 UX=AD*DX*U(I1)/UMAX UY=AD*DY*U(I2)/UMAX X2=X1+UX Y2=Y1+UY XX1=SQRT((UX*UX+UY*UY)/(DX*DX+DY*DY)) AR1=AR IF(IAR.EQ.1) AR1=AR*XX1*1.5 CALL GAROHD(X1,Y1,X2,Y2,AR1,CR) 12 CONTINUE Y1=Y1+DY DO 14 I=1,NX X1=X0+DX*FLOAT(NX+1-I) I1=NX+1-I+NX*(2*J-1) I2=I1+NXY R=SQRT(U(I1)*U(I1)+U(I2)*U(I2)) IF(R.LT.URMI) GO TO 14 UX=AD*DX*U(I1)/UMAX UY=AD*DY*U(I2)/UMAX X2=X1+UX Y2=Y1+UY XX1=SQRT((UX*UX+UY*UY)/(DX*DX+DY*DY)) AR1=AR IF(IAR.EQ.1) AR1=AR*XX1*1.5 CALL GAROHD(X1,Y1,X2,Y2,AR1,CR) 14 CONTINUE 10 CONTINUE RETURN C ENTRY GRAP2M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX,KP, & EXP) EXP=EXP NXY=NX*NY AD=0.5 DX=XL/FLOAT(NX+1) DY=YL/FLOAT(NY+1) DR=AMIN1(DX,DY) DRR=DR*1.0E-2 DRR=DR*2.0E-2 DRR=AMAX1(DRR,5.0E-3) UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 20 J=1,NY DO 20 I=1,NX Y=0.5*HY*FLOAT(2*J-NY-1) X=0.5*HX*FLOAT(2*I-NX-1+2*NXP) R=X*X+Y*Y IF(R.LT.AR2) GO TO 20 I1=I+NX*(J-1) X1=U(I1) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 20 CONTINUE IF(KP.EQ.0) GO TO 26 IF(KP.EQ.1) GO TO 22 IF(KP.EQ.-1) GO TO 24 UMIN=VMIN UMAX=VMAX GO TO 26 22 UMAX=VMAX GO TO 26 24 UMIN=VMIN 26 CONTINUE C DO 28 I=1,NXY IF(U(I).GT.UMAX) U(I)=UMAX IF(U(I).LT.UMIN) U(I)=UMIN 28 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,KP,UMIN,UMAX) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) Y=Y0+0.5*YL CALL PLOT(X,Y,3) X=X0+XL CALL PLOT(X,Y,2) X=X0+XEP*XL Y=Y0 CALL PLOT(X,Y,3) Y=Y0+YL CALL PLOT(X,Y,2) C DI=0.0 NY1=NY/2 DO 30 J=1,NY1 Y1=Y0+DY*FLOAT(2*J-1) DO 32 I=1,NX X1=X0+DX*FLOAT(I) I1=I+NX*(2*J-2) RAD=AD*DR*(U(I1)-UMIN)/(UMAX-UMIN) IF(RAD.LT.DRR) GO TO 32 X2=X1 Y2=Y1 CALL GCIRCL(X2,Y2,RAD) 32 CONTINUE C Y1=Y1+DY DO 34 I=1,NX X1=X0+DX*FLOAT(NX+1-I) I1=NX+1-I+NX*(2*J-1) RAD=AD*DR*(U(I1)-UMIN)/(UMAX-UMIN) IF(RAD.LT.DRR) GO TO 34 X2=X1 Y2=Y1 CALL GCIRCL(X2,Y2,RAD) 34 CONTINUE 30 CONTINUE RETURN C ENTRY GRAP3M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX, & V0,KP,INL,XEP) NXY=NX*NY AD=0.5 DX=XL/FLOAT(NX+1) DY=YL/FLOAT(NY+1) RR=DY/DX DR=AMIN1(DX,DY) DRR=DR*1.0E-2 DRR=DR*2.0E-2 DRR=AMAX1(DRR,5.0E-3) UMAX=U(1)+1.0E-10 UMIN=U(1)-1.0E-10 DO 40 J=1,NY DO 40 I=1,NX Y=0.5*HY*FLOAT(2*J-NY-1) X=0.5*HX*FLOAT(2*I-NX-1+2*NXP) R=X*X+Y*Y IF(R.LT.AR2) GO TO 40 I1=I+NX*(J-1) X1=U(I1) UMAX=AMAX1(UMAX,X1) UMIN=AMIN1(UMIN,X1) 40 CONTINUE IF(KP.EQ.0) V0=0.5*(UMIN+UMAX) IF(KP.EQ.0) GO TO 46 IF(KP.EQ.1) GO TO 42 IF(KP.EQ.-1) GO TO 44 UMIN=VMIN UMAX=VMAX GO TO 46 42 UMAX=VMAX GO TO 46 44 UMIN=VMIN 46 CONTINUE C DO 48 I=1,NXY IF(U(I).GT.UMAX) U(I)=UMAX IF(U(I).LT.UMIN) U(I)=UMIN 48 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,KP,UMIN,UMAX) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) Y=Y0+0.5*YL CALL PLOT(X,Y,3) X=X0+XL CALL PLOT(X,Y,2) X=X0+XEP*XL Y=Y0 CALL PLOT(X,Y,3) Y=Y0+YL CALL PLOT(X,Y,2) C DI=0.0 NY1=NY/2 X=ABS(UMAX-V0) Y=ABS(V0-UMIN) C DO 50 J=1,NY1 Y1=Y0+DY*FLOAT(2*J-1) DO 52 I=1,NX X1=X0+DX*FLOAT(I) I1=I+NX*(2*J-2) IF(U(I1).LT.V0) GO TO 54 RAD=AD*DR*(U(I1)-V0)/X IF(RAD.LT.DRR) GO TO 52 X2=X1 Y2=Y1 CALL GREC1(X2,Y2,RAD,RR,1) GO TO 52 54 CONTINUE RAD=AD*DR*(V0-U(I1))/Y IF(RAD.LT.DRR) GO TO 52 X2=X1 Y2=Y1 CALL GREC1(X2,Y2,RAD,RR,INL) 52 CONTINUE C Y1=Y1+DY DO 56 I=1,NX X1=X0+DX*FLOAT(NX+1-I) I1=NX+1-I+NX*(2*J-1) IF(U(I1).LT.V0) GO TO 58 RAD=AD*DR*(U(I1)-V0)/X IF(RAD.LT.DRR) GO TO 56 X2=X1 Y2=Y1 CALL GREC1(X2,Y2,RAD,RR,1) GO TO 56 58 CONTINUE RAD=AD*DR*(V0-U(I1))/Y IF(RAD.LT.DRR) GO TO 56 X2=X1 Y2=Y1 CALL GREC1(X2,Y2,RAD,RR,INL) 56 CONTINUE 50 CONTINUE RETURN END SUBROUTINE GBATU(X,Y,A,CR) ACR=A*CR X1=X-A Y1=Y-ACR X2=X+A Y2=Y+ACR CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) CALL PLOT(X2,Y1,3) CALL PLOT(X1,Y2,2) RETURN END SUBROUTINE GKOME(X,Y,A,CR) ACR=A*CR X1=X-A Y1=Y-ACR X2=X+A Y2=Y+ACR CALL PLOT(X1,Y1,3) CALL PLOT(X2,Y2,2) CALL PLOT(X2,Y1,3) CALL PLOT(X1,Y2,2) CALL PLOT(X,Y1,3) CALL PLOT(X,Y2,2) CALL PLOT(X2,Y,3) CALL PLOT(X1,Y,2) RETURN END SUBROUTINE ODATA(X,Y,LY,LM,LD,LH) C X,Y ------- DATE NO HYOGI ITI XM=0.5 YM=0.125 XX=999.0 X1=X Z=FLOAT(LY) CALL SYMBOL(X1,Y,YM,5HLAST=,0.0,5) CALL NUMBER(XX,XX,YM,Z,0.0,-1) Z=FLOAT(LM) CALL SYMBOL(XX,XX,YM,5H II=,0.0,5) CALL NUMBER(XX,XX,YM,Z,0.0,-1) Z=FLOAT(LD) CALL SYMBOL(XX,XX,YM,6H NXP=,0.0,6) CALL NUMBER(XX,XX,YM,Z,0.0,-1) Z=FLOAT(LH) CALL SYMBOL(XX,XX,YM,5H NR=,0.0,5) CALL NUMBER(XX,XX,YM,Z,0.0,-1) RETURN END SUBROUTINE OCHARA1(X,Y,IS,IO,IC,I,IA,UMIN,UMAX) IS=IS IO=IO IC=IC XX=999.0 YM=0.085 RIA=FLOAT(IA) IF(I.EQ.0) GO TO 20 CALL SYMBOL(X,Y,YM,3HII=,0.0,3) CALL NUMBER(XX,XX,YM,RIA,0.0,-1) CALL SYMBOL(XX,XX,YM,6H MIN=,0.0,6) CALL NUMBER(XX,XX,YM,UMIN,0.0,8) CALL SYMBOL(XX,XX,YM,6H MAX=,0.0,6) CALL NUMBER(XX,XX,YM,UMAX,0.0,8) GO TO 30 20 CONTINUE CALL SYMBOL(X,Y,YM,4HMIN=,0.0,4) CALL NUMBER(XX,XX,YM,UMIN,0.0,8) CALL SYMBOL(XX,XX,YM,6H MAX=,0.0,6) CALL NUMBER(XX,XX,YM,UMAX,0.0,8) 30 CONTINUE RETURN END FUNCTION CVMGP(X,Y,Z) CVMGP=Y IF(Z.GE.0.0) CVMGP=X RETURN END FUNCTION CVMGM(X,Y,Z) CVMGM=Y IF(Z.LT.0.0) CVMGM=X RETURN END SUBROUTINE GRAP8G(IA,LB,LL,NX,NXX,X0,Y0,XL,YL,U,VMIN,VMAX,KP) DIMENSION U(1) C CALL NEWPE1(2) IPD=0 I=LB-1 I=MOD(I,8) IF(I.EQ.0) IPD=1 XX=999.0 XM=0.075 YM=0.075 C NX1=NX+1 NX2=NX+2 N1=NXX+2 C INX=NXX/NX L1=LB+LL-1 C DO 8 J=LB,L1 C DO 8 I=1,NX C I1=I+N1*(J-1) C I2=I*INX+N1*(J-1) C U(I1)=U(I2) C 8 CONTINUE DX=XL/FLOAT(NX1) I1=1+N1*(LB-1) UMAX=U(I1)+1.0E-6 UMIN=U(I1)-1.0E-6 DO 10 J=LB,L1 DO 10 I=1,NX I1=I+N1*(J-1) UMAX=AMAX1(UMAX,U(I1)) UMIN=AMIN1(UMIN,U(I1)) 10 CONTINUE IF(KP.EQ.0) GO TO 16 IF(KP.EQ.1) GO TO 12 IF(KP.EQ.-1) GO TO 14 IF(KP.EQ.5) GO TO 15 UMIN=VMIN UMAX=VMAX GO TO 16 12 UMAX=VMAX GO TO 16 14 UMIN=VMIN GO TO 16 15 UMIN=-UMIN UMAX=AMAX1(UMAX,UMIN) UMIN=-UMAX 16 CONTINUE C DO 18 J=LB,L1 DO 18 I=1,NX I1=I+N1*(J-1) IF(U(I1).GT.UMAX) U(I1)=UMAX IF(U(I1).LT.UMIN) U(I1)=UMIN 18 CONTINUE C X1=X0 Y1=Y0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,1,IA,UMIN,UMAX) C CALL PLOT(X0,Y0,3) X=X0+XL Y=Y0 CALL PLOT(X,Y,2) Y=Y0+YL CALL PLOT(X,Y,2) X=X0 CALL PLOT(X,Y,2) Y=Y0 CALL PLOT(X,Y,2) C DO 20 JJ=LB,L1 I=4-MOD(JJ-LB,3)-IPD CALL NEWPE1(I) X=X0 I1=1+N1*(JJ-1) Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,3) DO 30 I=1,NX I1=I+N1*(JJ-1) X=X+DX Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,2) 30 CONTINUE I1=NX+N1*(JJ-1) X=X+DX Y=(U(I1)-UMIN)/(UMAX-UMIN) Y=Y0+YL*Y CALL PLOT(X,Y,2) 20 CONTINUE RETURN END c C FILE NAME /POST/SUBPS2 C FILE NAME / PSSC / SUBPS2 K. OHTA 1992/02/17 c file name /test3/subc.f subroutine xyopen write(6,10) 10 format( "%!" ) write(6,20) 20 format( "%%Page: 1 1" ) write(6,*) "/m { moveto } def" write(6,*) "/l { lineto } def" write(6,*) "/s { show } def" write(6,*) "/aws { awidthshow } def" write(6,*) "/sl { setlinewidth } def" write(6,*) "/sn { stroke newpath } def" write(6,*) "/tr { /Times-Roman findfont } def" write(6,*) "/sb { /Symbol findfont } def" write(6,*) "/sf { scalefont } def" write(6,*) "/se { setfont } def" write(6,*) "/ro { rotate } def" write(6,*) "/sd { setdash } def" write(6,*) "/st { stroke } def" write(6,*) "/sp { stroke showpage } def" write(6,*) "/tl { translate } def" write(6,*) "/sc { scale } def" write(6,*) "/np { newpath } def" write(6,*) "/cp { closepath } def" write(6,*) "/gs { gsave } def" write(6,*) "/gr { grestore } def" write(6,*) "/grc { grestore 0.0 0.0 0.0 setrgbcolor } def" write(6,*) "/sg { setgray } def" write(6,*) "/chs { sethsbcolor } def" write(6,*) "/crg { setrgbcolor } def" write(6,*) "/f { fill } def" write(6,*) "/cup { currentpoint } def" write(6,*) "/ftr /Times-Roman findfont 20 scalefont def" write(6,*) "/fnr /Helvetica findfont 20 scalefont def" write(6,*) "/fni /Helvetica-Oblique findfont 20 scalefont def" write(6,*) "/fnb /Helvetica-Bold findfont 30 scalefont def" write(6,*) "/shtr {moveto ftr setfont show} def" write(6,*) "/shnr {moveto fnr setfont show} def" write(6,*) "/shni {moveto fni setfont show} def" write(6,*) "/shnb {moveto fnb setfont show} def" write(6,*) "/c {setrgbcolor 1 0 translate newpath 0 0 moveto & 0 1 lineto 1 1 lineto 1 0 lineto closepath & fill} def" write(6,*) "/ch {sethsbcolor 1 0 translate newpath 0 0 moveto & 0 1 lineto 1 1 lineto 1 0 lineto closepath & fill} def" return end subroutine plots(idci,ity,isz,ino) c idci=1 dot, idci=2 cm, idci=3 in, ity=1 tate, ity=2 yoko, c isz=1 a4, isz=2 b4, ino=number common /cdci/ldci ldci=idci if(ino.ge.2) write(6,10) ino,ino 10 format("%%Page:",1x,i2,1x,i2) c c write(6,*) " 0.90 0.90 sc" c if(idci.eq.1) then if(ity.eq.1.and.(isz.eq.1.or.isz.eq.2)) then write(6,*) " 1.00 1.00 sc" write(6,*) " 28.0 33.0 tl" else if(ity.eq.2.and.isz.eq.1) then write(6,*) " 1.00 1.00 sc" write(6,*) " 90.0 ro" write(6,*) " 35.0 -565.0 tl" else if(ity.eq.2.and.isz.eq.2) then write(6,*) " 1.00 1.00 sc" write(6,*) " 90.0 ro" write(6,*) " 35.0 -696.0 tl" end if else if(idci.eq.2) then if(ity.eq.1.and.(isz.eq.1.or.isz.eq.2)) then write(6,*) " 28.38 28.17 sc" write(6,*) " 1.00 1.20 tl" else if(ity.eq.2.and.isz.eq.1) then write(6,*) " 90.0 ro" write(6,*) " 28.17 28.38 sc" write(6,*) " 1.15 -19.90 tl" else if(ity.eq.2.and.isz.eq.2) then write(6,*) " 90.0 ro" write(6,*) " 28.17 28.38 sc" write(6,*) " 1.15 -24.60 tl" end if else if(idci.eq.3) then if(ity.eq.1.and.(isz.eq.1.or.isz.eq.2)) then write(6,*) " 72.08 71.55 sc" write(6,*) " 0.39 0.47 tl" else if(ity.eq.2.and.isz.eq.1) then write(6,*) " 90.0 ro" write(6,*) " 71.55 72.08 sc" write(6,*) " 0.49 -7.88 tl" else if(ity.eq.2.and.isz.eq.2) then write(6,*) " 90.0 ro" write(6,*) " 71.55 72.08 sc" write(6,*) " 0.49 -9.66 tl" end if end if return end subroutine plote write(6,*) "sp" return end subroutine linee write(6,*) "st" return end subroutine dashe write(6,*) "sn [ ] 0 sd" return end subroutine plot(x1,y1,im) common /cdci/ldci if(ldci.eq.1) then if(im.eq.3) then write(6,10) x1,y1 10 format(f6.1,1x,f6.1," m") else if(im.eq.2) then write(6,20) x1,y1 20 format(f6.1,1x,f6.1," l") end if else if(ldci.ge.2) then if(im.eq.3) then write(6,30) x1,y1 30 format(f6.2,1x,f6.2," m") else if(im.eq.2) then write(6,40) x1,y1 40 format(f6.2,1x,f6.2," l") end if end if return end subroutine factor(fct) write(6,10) fct,fct 10 format(f6.3,1x,f6.3," sc") return end subroutine newpen(ipn) common /cdci/ldci if(ldci.eq.1) pn=0.350 if(ldci.eq.2) pn=0.012 if(ldci.eq.3) pn=0.005 ip=mod(ipn-1,10)+1 is=mod((ipn/10)-1,10)+1 ip=max0(ip,1) fp=float(ip)*pn write(6,30) fp 30 format("sn",1x,f6.3," sl") if(ipn.gt.10) then if(ldci.eq.1) then if(is.eq.1) write(6,*) "[50.0 0.00] 0 sd" if(is.eq.2) write(6,*) "[1.50 1.50] 0 sd" if(is.eq.3) write(6,*) "[3.00 3.00] 0 sd" if(is.eq.4) write(6,*) "[6.00 6.00] 0 sd" if(is.eq.5) write(6,*) "[9.00 9.00] 0 sd" if(is.eq.6) write(6,*) "[12.0 12.0] 0 sd" if(is.eq.7) write(6,*) "[6.00 3.00] 0 sd" if(is.eq.8) write(6,*) "[12.0 3.00] 0 sd" if(is.eq.9) write(6,*) "[18.0 3.00] 0 sd" if(is.eq.10) write(6,*) "[24.0 3.00] 0 sd" else if(ldci.eq.2) then if(is.eq.1) write(6,*) "[1.00 0.00] 0 sd" if(is.eq.2) write(6,*) "[0.05 0.05] 0 sd" if(is.eq.3) write(6,*) "[0.10 0.10] 0 sd" if(is.eq.4) write(6,*) "[0.20 0.20] 0 sd" if(is.eq.5) write(6,*) "[0.30 0.30] 0 sd" if(is.eq.6) write(6,*) "[0.40 0.40] 0 sd" if(is.eq.7) write(6,*) "[0.20 0.10] 0 sd" if(is.eq.8) write(6,*) "[0.40 0.10] 0 sd" if(is.eq.9) write(6,*) "[0.60 0.10] 0 sd" if(is.eq.10) write(6,*) "[0.80 0.10] 0 sd" else if(ldci.eq.3) then if(is.eq.1) write(6,*) "[1.00 0.00] 0 sd" if(is.eq.2) write(6,*) "[0.02 0.02] 0 sd" if(is.eq.3) write(6,*) "[0.04 0.04] 0 sd" if(is.eq.4) write(6,*) "[0.08 0.08] 0 sd" if(is.eq.5) write(6,*) "[0.12 0.12] 0 sd" if(is.eq.6) write(6,*) "[0.16 0.16] 0 sd" if(is.eq.7) write(6,*) "[0.08 0.04] 0 sd" if(is.eq.8) write(6,*) "[0.16 0.04] 0 sd" if(is.eq.9) write(6,*) "[0.24 0.04] 0 sd" if(is.eq.10) write(6,*) "[0.32 0.04] 0 sd" end if end if return end subroutine transl(x,y) write(6,22) x,y 22 format(f7.2,1x,f7.2," tl") return end subroutine symbol(x,y,h,isymb,ang,n) character isymb*80,ica*80,ich(80)*1 equivalence (ica,ich) ica=isymb write(6,*) "gs" h0=h*1.50 write(6,10) h0 10 format("tr ",f6.2," sf se") if(x.gt.998.0.or.y.gt.998.0) then write(6,*) "cup m" else write(6,20) x,y 20 format(f6.2,1x,f6.2," m") end if write(6,30) ang 30 format(f6.1," ro") h1=h*0.06 h2=h*0.02 write(6,40) h1,h2 40 format(f6.2," 0.0 8#040 ",f6.2," 0.0") c write(6,*) "(",(ich(i),i=1,n),") s" write(6,*) "(",(ich(i),i=1,n),") aws" write(6,*) "gr" return end subroutine symblc(x,y,h,isymb,ang,n) character isymb*80,ica*80,ich(80)*1 equivalence (ica,ich) ica=isymb write(6,*) "gs" h0=h*1.50 write(6,10) h0 10 format("sb ",f6.2," sf se") if(x.gt.998.0.or.y.gt.998.0) then write(6,*) "cup m" else write(6,20) x,y 20 format(f6.2,1x,f6.2," m") end if write(6,30) ang 30 format(f6.1," ro") h1=h*0.06 h2=h*0.02 write(6,40) h1,h2 40 format(f6.2," 0.0 8#040 ",f6.2," 0.0") c write(6,*) "(",(ich(i),i=1,n),") s" write(6,*) "(",(ich(i),i=1,n),") aws" write(6,*) "gr" return end subroutine circ1(x,y,r) x1=x+r write(6,10) x1,y 10 format(f6.2,1x,f6.2," m") write(6,20) x,y,r 20 format(3(1x,f6.2)," 0 360 arc") return end subroutine rect(x,y,h,w,ang,ipn) ar=ang*3.1415926/180.0 as=sin(ar) ac=cos(ar) xw=w*ac yw=w*as xh=h*as yh=h*ac x1=x y1=y x2=x-xh y2=y+yh x3=x+xw-xh y3=y+yw+yh x4=x+xw y4=y+yw call plot(x1,y1,ipn) call plot(x2,y2,2) call plot(x3,y3,2) call plot(x4,y4,2) call plot(x1,y1,2) return end subroutine grab write(6,*) "gs" return end subroutine grae write(6,*) "gr" return end subroutine newhsb(h,s,b) ah=amax1(h,0.0) as=amax1(s,0.0) ab=amax1(b,0.0) ah=amin1(h,1.0) as=amin1(s,1.0) ab=amin1(b,1.0) write(6,10) ah,as,ab 10 format(3(1x,f4.2)," chs") return end subroutine newrgb(r,g,b) ar=amax1(r,0.0) ag=amax1(g,0.0) ab=amax1(b,0.0) ar=amin1(r,1.0) ag=amin1(g,1.0) ab=amin1(b,1.0) write(6,10) ar,ag,ab 10 format(3(1x,f4.2)," crg") return end subroutine newgry(ig) is=max0(ig,1) is=min0(ig,10) fg=float(is-1)*0.1 write(6,10) fg 10 format(f4.2," sg") return end subroutine paint write(6,*) "f" return end subroutine number(x,y,h,fnp,ang,n) character fmt1*6,fmt2*7,fmt3*4,fmt4*5,fna*80 fnp1=fnp fnpa=abs(fnp1)+1.0e-30 keta=1 if(fnpa.gt.1.0) keta=ifix(alog10(fnpa)*1.00001)+1 if(fnpa.lt.0.01) keta=-(ifix(alog10(fnpa)*1.00001))-1 if(fnpa.lt.1.0e-9) keta=1 if(abs(fnp).lt.1.0e-9) fnp1=1.0e-9 minus=0 if(fnp1.lt.0.0) minus=1 if(n.lt.0) go to 1 n1=keta+n+1+minus if(n1.le.0) n1=1 if(n1.ge.10) go to 2 write(fmt1,100) n1,n 100 format('(f',i1,'.',i1,')') write(fna,fmt1) fnp1 call symbol(x,y,h,fna,ang,n1) 2 if(n1.lt.10) go to 1 write(fmt2,110) n1,n 110 format('(f',i2,'.',i1,')') write(fna,fmt2) fnp1 call symbol(x,y,h,fna,ang,n1) c 1 if(n.ge.0) go to 3 if(fnp1.gt.0.0) ifnp=ifix(fnp1*(10.0**(n+1))+0.5) if(fnp1.lt.0.0) ifnp=ifix(fnp1*(10.0**(n+1))-0.5) n2=keta+n+1+minus if(n2.le.0) n2=1 if(n2.ge.10) go to 4 write(fmt3,200) n2 200 format('(i',i1,')') write(fna,fmt3) ifnp call symbol(x,y,h,fna,ang,n2) 4 if(n2.lt.10) go to 3 write(fmt4,210) n2 210 format('(i',i2,')') write(fna,fmt4) ifnp call symbol(x,y,h,fna,ang,n2) 3 return end subroutine color1(nx,ny,xb,yb,xl,yl,u,ico) dimension u(1) dx=xl/float(nx-1) dy=yl/float(ny-1) c write(6,*) "gs" write(6,20) dx,dy 20 format(f7.4,1x,f7.4," sc") c ya=yb/dy-1.5 xa=xb/dx-1.5+float(nx) write(6,22) xa,ya c do 100 j=1,ny y=1.0 x=-float(nx) xx=x yy=y write(6,22) xx,yy 22 format(f7.2,1x,f7.2," tl") do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 r=a0-(a0-0.5)*x1 g=a0-(a0-0.5)*x1 b=a0-(a0-0.5)*x1 go to 40 32 continue r=a0-a0*x1 g=1.0 b=a0-a0*x1 go to 40 33 continue r=1.0 g=a0-a0*x1 b=a0-a0*x1 go to 40 34 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) g=a0+3.0-y1 b=a0+1.0-y1 if(y1.ge.a0+2.0) b=y1-(a0+3.0) go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) write(6,24) r,g,b 24 format(f4.2,1x,f4.2,1x,f4.2," c") c* 24 format(f4.2,1x,f4.2,1x,f4.2," ch") 100 continue c x1=-float(nx-1) y1=-float(ny-1) write(6,22) x1,y1 c x1=0.5 y1=0.5 x2=x1+float(nx-1) y2=y1+float(ny-1) call newrgb(1.0,0.0,0.0) call plot(x1,y1,3) call plot(x1,y2,2) call plot(x2,y2,2) call plot(x2,y1,2) call plot(x1,y1,2) write(6,*) "st" write(6,*) "grc" return end subroutine color2(nx,ny,xb,yb,xl,yl,u,vmin,vmax,ico) dimension u(1) dx=xl/float(nx-1) dy=yl/float(ny-1) c do 10 j=1,ny do 10 i=1,nx i1=i+nx*(j-1) u(i1)=(u(i1)-vmin)/(vmax-vmin) 10 continue c write(6,*) "gs" write(6,20) dx,dy 20 format(f7.4,1x,f7.4," sc") c ya=yb/dy-1.5 xa=xb/dx-1.5+float(nx) write(6,22) xa,ya c do 100 j=1,ny y=1.0 x=-float(nx) xx=x yy=y write(6,22) xx,yy 22 format(f7.2,1x,f7.2," tl") do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 c x1=2.0*x1 r=x1 g=2.0-x1 b=0.0 go to 40 32 continue r=x1 g=x1 b=x1 go to 40 33 continue x1=7.0*x1 b=x1 if(x1.ge.2.0) b=3.0-x1 if(x1.ge.5.0) b=0.5*(x1-5.0) g=x1-1.0 if(x1.ge.4.0) g=5.0-x1 if(x1.ge.5.0) g=0.5*(x1-5.0) r=x1-3.0 go to 40 34 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) write(6,24) r,g,b 24 format(f4.2,1x,f4.2,1x,f4.2," c") c* 24 format(f4.2,1x,f4.2,1x,f4.2," ch") 100 continue c x1=-float(nx-1) y1=-float(ny-1) write(6,22) x1,y1 c x1=0.5 y1=0.5 x2=x1+float(nx-1) y2=y1+float(ny-1) x12=0.5*(x1+x2) y12=0.5*(y1+y2) x13=x1+(61.0/321.0)*(x2-x1) c call newrgb(1.0,0.0,0.0) call newrgb(0.0,0.0,1.0) call plot(x1,y1,3) call plot(x1,y2,2) call plot(x2,y2,2) call plot(x2,y1,2) call plot(x1,y1,2) call plot(x1,y12,3) call plot(x2,y12,2) if(nx.eq.ny) then call plot(x12,y1,3) call plot(x12,y2,2) else call plot(x13,y1,3) call plot(x13,y2,2) endif write(6,*) "st" write(6,*) "grc" return end subroutine imagc1(nx,ny,xb,yb,xl,yl,vmin,vmax,u,ico) parameter (nc=60) dimension u(1),ird(nc),igd(nc),ibd(nc) c write(6,*) "gs" write(6,22) xb,yb 22 format(f7.2,1x,f7.2," tl") write(6,20) xl,yl 20 format(f7.4,1x,f7.4," sc") c do 100 k=1,3 do 100 j=1,ny do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 if(ico.eq.7) go to 37 r=a0-(a0-0.5)*x1 g=a0-(a0-0.5)*x1 b=a0-(a0-0.5)*x1 go to 40 32 continue r=a0-a0*x1 g=1.0 b=a0-a0*x1 go to 40 33 continue r=1.0 g=a0-a0*x1 b=a0-a0*x1 go to 40 34 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) g=a0+3.0-y1 b=a0+1.0-y1 if(y1.ge.a0+2.0) b=y1-(a0+3.0) go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 37 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c i2=mod(i1-1,nc)+1 ird(i2)=r*15.0 igd(i2)=g*15.0 ibd(i2)=b*15.0 c if(k.eq.1) then if(j.eq.1.and.i.eq.1) write(6,51) 51 format("/Red <") if(mod(i1,nc).eq.0) write(6,52) ird 52 format(60z1) if(j.eq.ny.and.i.eq.nx) write(6,53) ird 53 format(60z1,/,"> def") else if(k.eq.2) then if(j.eq.1.and.i.eq.1) write(6,55) 55 format("/Green <") if(mod(i1,nc).eq.0) write(6,52) igd if(j.eq.ny.and.i.eq.nx) write(6,53) igd else if(k.eq.3) then if(j.eq.1.and.i.eq.1) write(6,57) 57 format("/Blue <") if(mod(i1,nc).eq.0) write(6,52) ibd if(j.eq.ny.and.i.eq.nx) write(6,53) ibd endif 100 continue c write(6,61) nx,ny 61 format(i3,1x,i3," 4") c write(6,62) nx,-ny,ny c 62 format("[",i3," 0 0 ",i4," 0 ",i3,"]") write(6,62) nx,ny 62 format("[",i3," 0 0 ",i3," 0 0]") write(6,63) 63 format("{Red} {Green} {Blue}",/,"true 3 colorimage") c x1=0.0 y1=0.0 x2=1.0 y2=1.0 c call newpen(1) c call newrgb(1.0,0.0,0.0) c call plot(x1,y1,3) c call plot(x1,y2,2) c call plot(x2,y2,2) c call plot(x2,y1,2) c call plot(x1,y1,2) write(6,*) "st" write(6,*) "grc" return end subroutine image1(nx,ny,xb,yb,xl,yl,u,vmin,vmax,ico) parameter (nc=60) dimension u(1),iwd(nc) c do 10 j=1,ny do 10 i=1,nx i1=i+nx*(j-1) u(i1)=(u(i1)-vmin)/(vmax-vmin) 10 continue c write(6,*) "gs" write(6,22) xb,yb 22 format(f7.2,1x,f7.2," tl") write(6,20) xl,yl 20 format(f7.4,1x,f7.4," sc") c a0=1.05 b0=0.05 c write(6,*) "<" do 100 j=1,ny do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c i2=mod(i1-1,nc)+1 if(ic.lt.0) w=a0-(a0-b0)*x1 if(ic.gt.0) w=b0+(a0-b0)*x1 w=amax1(w,0.0) w=amin1(w,1.0) iwd(i2)=w*15.0 if(mod(i1,nc).eq.0) write(6,52) iwd if(j.eq.ny.and.i.eq.nx) write(6,52) iwd 52 format(60z1) c 100 continue c write(6,*) ">" write(6,61) nx,ny 61 format(i3,1x,i3," 4") c write(6,62) nx,-ny,ny c 62 format("[",i3," 0 0 ",i4," 0 ",i3,"]") write(6,62) nx,ny 62 format("[",i3," 0 0 ",i3," 0 0]") write(6,63) 63 format("{ } image") write(6,*) "grc" c c x1=-float(nx-1) y1=-float(ny-1) c write(6,22) x1,y1 c x1=xb y1=yb x2=xb+xl y2=yb+yl x12=0.5*(x1+x2) y12=0.5*(y1+y2) x13=x1+(61.0/321.0)*(x2-x1) c call newrgb(1.0,0.0,0.0) c call newrgb(0.0,0.0,1.0) call newrgb(0.0,0.0,0.0) call plot(x1,y1,3) call plot(x1,y2,2) call plot(x2,y2,2) call plot(x2,y1,2) call plot(x1,y1,2) call plot(x1,y12,3) call plot(x2,y12,2) if(nx.eq.ny) then call plot(x12,y1,3) call plot(x12,y2,2) else call plot(x13,y1,3) call plot(x13,y2,2) endif write(6,*) "st" c call newrgb(0.0,0.0,1.0) call newrgb(0.0,0.0,0.0) c write(6,*) "grc" return end subroutine imagc2(nx,ny,xb,yb,xl,yl,u,vmin,vmax,ico) parameter (nc=60) dimension u(1),ird(nc),igd(nc),ibd(nc) c do 10 j=1,ny do 10 i=1,nx i1=i+nx*(j-1) u(i1)=(u(i1)-vmin)/(vmax-vmin) 10 continue c write(6,*) "gs" write(6,22) xb,yb 22 format(f7.2,1x,f7.2," tl") write(6,20) xl,yl 20 format(f7.4,1x,f7.4," sc") c do 100 k=1,3 do 100 j=1,ny do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 c x1=2.0*x1 r=x1 g=2.0-x1 b=0.0 go to 40 32 continue r=x1 g=x1 b=x1 r=4.0*(x1-0.5) b=4.0*(0.5-x1) if(x1.le.0.25) g=4.0*(0.25-x1) if(x1.gt.0.25) g=4.0*(x1-0.75) go to 40 c r=x1 c g=x1 c b=x1 c go to 40 33 continue x1=7.0*x1 b=x1 if(x1.ge.2.0) b=3.0-x1 if(x1.ge.5.0) b=0.5*(x1-5.0) g=x1-1.0 if(x1.ge.4.0) g=5.0-x1 if(x1.ge.5.0) g=0.5*(x1-5.0) r=x1-3.0 go to 40 34 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c i2=mod(i1-1,nc)+1 ird(i2)=r*15.0 igd(i2)=g*15.0 ibd(i2)=b*15.0 c if(k.eq.1) then if(j.eq.1.and.i.eq.1) write(6,51) 51 format("/Red <") if(mod(i1,nc).eq.0) write(6,52) ird 52 format(60z1) if(j.eq.ny.and.i.eq.nx) write(6,53) ird 53 format(60z1,/,"> def") else if(k.eq.2) then if(j.eq.1.and.i.eq.1) write(6,55) 55 format("/Green <") if(mod(i1,nc).eq.0) write(6,52) igd if(j.eq.ny.and.i.eq.nx) write(6,53) igd else if(k.eq.3) then if(j.eq.1.and.i.eq.1) write(6,57) 57 format("/Blue <") if(mod(i1,nc).eq.0) write(6,52) ibd if(j.eq.ny.and.i.eq.nx) write(6,53) ibd endif 100 continue c write(6,61) nx,ny 61 format(i3,1x,i3," 4") c write(6,62) nx,-ny,ny c 62 format("[",i3," 0 0 ",i4," 0 ",i3,"]") write(6,62) nx,ny 62 format("[",i3," 0 0 ",i3," 0 0]") write(6,63) 63 format("{Red} {Green} {Blue}",/,"true 3 colorimage") c x1=0.0 y1=0.0 x2=1.0 y2=1.0 c call newpen(1) c call newrgb(1.0,0.0,0.0) c call plot(x1,y1,3) c call plot(x1,y2,2) c call plot(x2,y2,2) c call plot(x2,y1,2) c call plot(x1,y1,2) c write(6,*) "st" write(6,*) "grc" c x1=-float(nx-1) y1=-float(ny-1) c write(6,22) x1,y1 c x1=xb y1=yb x2=xb+xl y2=yb+yl x12=0.5*(x1+x2) y12=0.5*(y1+y2) c x13=x1+(61.0/501.0)*(x2-x1) c x13=x1+(101.0/501.0)*(x2-x1) x13=x1+(61.0/181.0)*(x2-x1) c call newrgb(1.0,0.0,0.0) call newrgb(0.0,0.0,1.0) call plot(x1,y1,3) call plot(x1,y2,2) call plot(x2,y2,2) call plot(x2,y1,2) call plot(x1,y1,2) call plot(x1,y12,3) call plot(x2,y12,2) if(nx.eq.ny) then call plot(x12,y1,3) call plot(x12,y2,2) else call plot(x13,y1,3) call plot(x13,y2,2) endif write(6,*) "st" call newrgb(0.0,0.0,0.0) c write(6,*) "grc" c return end subroutine image2(nx,ny,xb,yb,xl,yl,u,vmin,vmax,ico) parameter (nc=30) dimension u(1),iwd(nc) c do 10 j=1,ny do 10 i=1,nx i1=i+nx*(j-1) u(i1)=(u(i1)-vmin)/(vmax-vmin) 10 continue c write(6,*) "gs" write(6,22) xb,yb 22 format(f7.2,1x,f7.2," tl") write(6,20) xl,yl 20 format(f7.4,1x,f7.4," sc") c a0=1.05 b0=0.05 c write(6,*) "<" do 100 j=1,ny do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c i2=mod(i1-1,nc)+1 if(ic.lt.0) w=a0-(a0-b0)*x1 if(ic.gt.0) w=b0+(a0-b0)*x1 w=amax1(w,0.0) w=amin1(w,1.0) iwd(i2)=w*255.0 if(mod(i1,nc).eq.0) write(6,52) iwd if(j.eq.ny.and.i.eq.nx) write(6,52) iwd 52 format(30z2.2) c 100 continue c write(6,*) ">" write(6,61) nx,ny 61 format(i3,1x,i3," 8") c write(6,62) nx,-ny,ny c 62 format("[",i3," 0 0 ",i4," 0 ",i3,"]") write(6,62) nx,ny 62 format("[",i3," 0 0 ",i3," 0 0]") write(6,63) 63 format("{ } image") write(6,*) "grc" c c x1=-float(nx-1) y1=-float(ny-1) c write(6,22) x1,y1 c x1=xb y1=yb x2=xb+xl y2=yb+yl x12=0.5*(x1+x2) y12=0.5*(y1+y2) x13=x1+(61.0/321.0)*(x2-x1) c x13=x1+(101.0/501.0)*(x2-x1) c call newrgb(1.0,0.0,0.0) c call newrgb(0.0,0.0,1.0) call newrgb(0.0,0.0,0.0) call plot(x1,y1,3) call plot(x1,y2,2) call plot(x2,y2,2) call plot(x2,y1,2) call plot(x1,y1,2) call plot(x1,y12,3) call plot(x2,y12,2) if(nx.eq.ny) then call plot(x12,y1,3) call plot(x12,y2,2) else call plot(x13,y1,3) call plot(x13,y2,2) endif write(6,*) "st" c call newrgb(0.0,0.0,1.0) call newrgb(0.0,0.0,0.0) c write(6,*) "grc" return end subroutine imagc8(nx,ny,xb,yb,xl,yl,u,vmin,vmax,ico) parameter (nc=30) dimension u(1),ird(nc),igd(nc),ibd(nc) c do 10 j=1,ny do 10 i=1,nx i1=i+nx*(j-1) u(i1)=(u(i1)-vmin)/(vmax-vmin) 10 continue c write(6,*) "gs" write(6,22) xb,yb 22 format(f7.2,1x,f7.2," tl") write(6,20) xl,yl 20 format(f7.4,1x,f7.4," sc") c do 100 k=1,3 do 100 j=1,ny do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 c x1=2.0*x1 r=x1 g=2.0-x1 b=0.0 go to 40 32 continue r=x1 g=x1 b=x1 go to 40 33 continue x1=7.0*x1 b=x1 if(x1.ge.2.0) b=3.0-x1 if(x1.ge.5.0) b=0.5*(x1-5.0) g=x1-1.0 if(x1.ge.4.0) g=5.0-x1 if(x1.ge.5.0) g=0.5*(x1-5.0) r=x1-3.0 go to 40 34 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c i2=mod(i1-1,nc)+1 ird(i2)=r*255.0 igd(i2)=g*255.0 ibd(i2)=b*255.0 c if(k.eq.1) then if(j.eq.1.and.i.eq.1) write(6,51) 51 format("/Red <") if(mod(i1,nc).eq.0) write(6,52) ird 52 format(30z2.2) if(j.eq.ny.and.i.eq.nx) write(6,53) ird 53 format(30z2.2,/,"> def") else if(k.eq.2) then if(j.eq.1.and.i.eq.1) write(6,55) 55 format("/Green <") if(mod(i1,nc).eq.0) write(6,52) igd if(j.eq.ny.and.i.eq.nx) write(6,53) igd else if(k.eq.3) then if(j.eq.1.and.i.eq.1) write(6,57) 57 format("/Blue <") if(mod(i1,nc).eq.0) write(6,52) ibd if(j.eq.ny.and.i.eq.nx) write(6,53) ibd endif 100 continue c write(6,61) nx,ny 61 format(i3,1x,i3," 8") c write(6,62) nx,-ny,ny c 62 format("[",i3," 0 0 ",i4," 0 ",i3,"]") write(6,62) nx,ny 62 format("[",i3," 0 0 ",i3," 0 0]") write(6,63) 63 format("{Red} {Green} {Blue}",/,"true 3 colorimage") c x1=0.0 y1=0.0 x2=1.0 y2=1.0 c call newpen(1) c call newrgb(1.0,0.0,0.0) c call plot(x1,y1,3) c call plot(x1,y2,2) c call plot(x2,y2,2) c call plot(x2,y1,2) c call plot(x1,y1,2) c write(6,*) "st" write(6,*) "grc" c x1=-float(nx-1) y1=-float(ny-1) c write(6,22) x1,y1 c x1=xb y1=yb x2=xb+xl y2=yb+yl x12=0.5*(x1+x2) y12=0.5*(y1+y2) x13=x1+(61.0/321.0)*(x2-x1) c call newrgb(1.0,0.0,0.0) call newrgb(0.0,0.0,1.0) call plot(x1,y1,3) call plot(x1,y2,2) call plot(x2,y2,2) call plot(x2,y1,2) call plot(x1,y1,2) call plot(x1,y12,3) call plot(x2,y12,2) if(nx.eq.ny) then call plot(x12,y1,3) call plot(x12,y2,2) else call plot(x13,y1,3) call plot(x13,y2,2) endif write(6,*) "st" call newrgb(0.0,0.0,0.0) c write(6,*) "grc" c return end