C FILE NAME / POSC / SUBPS10 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 c 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 xo=x yo=y zo=z if(yo.lt.0.0) then y=-y z=-z endif NX1=NX+1 NX2=NX+2 NY1=NY+1 NZ1=NZ+1 NY2=NY+2 NZ2=NZ+2 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) c XJ=Y/HY+0.5*FLOAT(ny2+1) xj=y/hy+1.5 XK=Z/HZ+1.5 c XK=Z/HZ+0.5*FLOAT(nz2+1) 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,3 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) P(I2+5)=P(I2) 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 if(yo.lt.0.0) then i1=1+L1 p(i1)=-p(i1) endif x=xo y=yo z=zo RETURN END SUBROUTINE QUANT2(X,Y,Z,HX,HY,HZ,MM,NX,NY,NZ,NXP,AA,F,P) DIMENSION F(1),P(1),AA(1),B(3) C xo=x yo=y zo=z if(yo.lt.0.0) then y=-y z=-z endif NX1=NX+1 NX2=NX+2 NY1=NY+1 ny2=ny+2 NZ1=NZ+1 nz2=nz+2 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 c XK=Z/HZ+1.5 xk=z/hz+0.5*float(nz2+1) 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,3 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) P(I2+5)=P(I2) 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 if(yo.lt.0.0) then i1=1+L1 p(i1)=-p(i1) endif x=xo y=yo z=zo RETURN END C FILE NAME /POSC/SUBPS11 SUBROUTINE GRAP7G(IA,LB,LL,NX,NXX,NXP,X0,Y0,XL,YL, & U,VMIN,VMAX,KP) DIMENSION U(1) C 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 NEWPEN(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,IPEN) C DIMENSION U(1) C CALL NEWPEN(1) 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) call linee C ULD=FLOAT(LANK+1)/(UMAX-UMIN) EP1=1.0E-8 CALL NEWPEN(IPEN) C DO 110 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 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 if(mod(i,50).eq.0) call linee 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) C DIMENSION U(1) PI=3.1415926 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=1 ILO=1 CALL NEWPEN(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 NEWPEN(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 C DO 110 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+(UMAX-UMIN)*ULD IF(U0.LT.X) GO TO 60 C CALL NEWPEN(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 NEWPEN(IHI) IF(U0.LT.U00) CALL NEWPEN(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 NEWPEN(IHI) IF(U0.LT.U00) CALL NEWPEN(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 NEWPEN(IHI) IF(U0.LT.U00) CALL NEWPEN(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 if(mod(i,50).eq.0) call linee C 100 CONTINUE call linee 110 CONTINUE CALL NEWPEN(IPEN) RETURN END SUBROUTINE GRAP3A(NR,NZ,MB,NX,NY,NXP,X0,Y0,XL,YL, & AR2,U,VMIN,VMAX,URMI0,KP,IAR,LAS1) C PARAMETER XM=3.4, YM=-0.6 DIMENSION U(1) C CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8 C XM=3.4 YM=-0.6 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) DX=XL/FLOAT(NR+1) DY=YL/FLOAT(NZ+1) AR=DX*HX*0.40 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 C 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 C NY1=NY/2+1 DO 10 J=1,NY 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 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) call linee 12 CONTINUE call linee 10 CONTINUE RETURN END SUBROUTINE AINTE1(IA,AA,F,P) PARAMETER(nline=2000,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt) 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) IPRO=IA(10) 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) ccc waku kaku tokoro 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 y=gy0 call plot(x,y,2) X=GX0 CALL PLOT(X,Y,2) 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 X=GX0+xdd1 Y=GY0+GYL*0.5 CALL PLOT(X,Y,2) Y=GY0 CALL PLOT(X,Y,2) c call linee ccc waku kakiowatta! C ccc boundary?? condition ?? 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 gymi=-gyma c gzmi=-gyma c GZMI GZMA ?? GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 C cc polar cap kara start suru basyo wo kimeru JB1=NY2/2+1 do 100 jj=1,2 DO 100 J=JB1,NY2 cc DO 100 J=1,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 c z ga + sika naiyouna.... Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=ARB*ARB-X1*X1-Y1*Y1 c zga 0.0 yori tiisaitoki 100 he IF(Z1.LE.0.0) GO TO 100 c IF(Z1.LE.0.0) GO TO 1030 if(jj.eq.1) Z1=SQRT(Z1) if(jj.eq.2) Z1=-SQRT(Z1) 1040 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) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) IF(IPRO.EQ.1) GX1=GX0+XDD1+AMP*X1 IF(IPRO.EQ.1) GY1=GY0+0.5*GYL+AMP*AMPY*Z1 IF(IPRO.EQ.2) GX1=GX0+XDD1+AMP*X1 IF(IPRO.EQ.2) GY1=GY0+0.5*GYL-AMP*AMPY*Y1 c CALL PLOT(GX1,GY1,3) lin=1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 c**** 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 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 c IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 IF(Z2.GT.GZMA) GO TO 200 IF(Z2.LT.GZMI) GO TO 210 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) R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 c IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 IF(Z3.GT.GZMA) GO TO 200 IF(Z3.LT.GZMI) GO TO 210 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) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) IF(IPRO.EQ.1) GX1=GX0+XDD1+AMP*X3 IF(IPRO.EQ.1) GY1=GY0+0.5*GYL+AMP*AMPY*Z3 IF(IPRO.EQ.2) GX1=GX0+XDD1+AMP*X3 IF(IPRO.EQ.2) GY1=GY0+0.5*GYL-AMP*AMPY*Y3 c CALL PLOT(GX1,GY1,2) lin=iii+1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,100).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c**** 200 continue if(lin.eq.1) go to 100 call newrgb(0.0,0.5,1.0) go to 300 210 continue if(lin.eq.1) go to 100 call newrgb(0.0,1.0,0.0) 300 continue lasl=lin lin=1 li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,3) do 40 ii=2,lasl lin=ii li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,2) if(mod(ii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 40 continue call linee 100 CONTINUE call newrgb(0.0,0.0,0.0) RETURN END SUBROUTINE AINTE2(IA,AA,F,P) PARAMETER(nline=2000,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt) 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) XLIM=AA(23) XLI2=AA(24) 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 C y=gy0 call plot(x,y,2) X=GX0 CALL PLOT(X,Y,2) 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 X=GX0+xdd1 Y=GY0+GYL*0.5 CALL PLOT(X,Y,2) Y=GY0 CALL PLOT(X,Y,2) c call linee 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=XXL/FLOAT(NX1) DY=YYL/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) gymi=-gyma AMPY=0.5*(GXMA-GXMI)*GYL/(GYMA*GXL) c GZMI=-0.1*HZ GZMI=-GYMA 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 jj=1,2 DO 100 J=JB1,NY2 C DO 100 J=1,NY1 c DO 100 J=1,NY2 DO 100 I=1,NX2 X1=0.5*DX*(2*I-NX2-1)+HX*FLOAT(NXP) Y1=0.5*DY*(2*J-NY2-1) C Y1=0.5*DY*(J) C y1=0.5*(2*J-1) Z1=0.0 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(M3) IF(X2.LT.0.0) AA(11)=-1.0 if(jj.eq.2) aa(11)=-aa(11) 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) c CALL PLOT(GX1,GY1,3) lin=1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 c**** 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 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 c XLIM=0.0 IF(X2.LT.XLIM.AND.Y2.LT.GYMI) GO TO 100 IF(X2.LT.XLIM.AND.Y2.GT.GYMA) GO TO 100 IF(X2.LT.XLIM.AND.Z2.LT.GZMI) GO TO 100 IF(X2.LT.XLIM.AND.Z2.GT.GZMA) GO TO 100 IF(X2.LT.XLI2.AND.Y2.LT.GYMI) GO TO 220 IF(X2.LT.XLI2.AND.Y2.GT.GYMA) GO TO 220 IF(X2.LT.XLI2.AND.Z2.LT.GZMI) GO TO 220 IF(X2.LT.XLI2.AND.Z2.GT.GZMA) GO TO 220 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 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) R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.XLIM.AND.Y3.LT.GYMI) GO TO 100 IF(X3.LT.XLIM.AND.Y3.GT.GYMA) GO TO 100 IF(X3.LT.XLIM.AND.Z3.LT.GZMI) GO TO 100 IF(X3.LT.XLIM.AND.Z3.GT.GZMA) GO TO 100 IF(X3.LT.XLI2.AND.Y3.LT.GYMI) GO TO 220 IF(X3.LT.XLI2.AND.Y3.GT.GYMA) GO TO 220 IF(X3.LT.XLI2.AND.Z3.LT.GZMI) GO TO 220 IF(X3.LT.XLI2.AND.Z3.GT.GZMA) GO TO 220 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 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) c CALL PLOT(GX1,GY1,2) lin=iii+1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,100).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c go to 100 c**** 200 continue if(lin.eq.1) go to 100 call newrgb(1.0,0.0,0.0) go to 300 220 continue if(lin.eq.1) go to 100 call newrgb(1.0,0.8,0.0) go to 300 210 continue if(lin.eq.1) go to 100 call newrgb(0.0,1.0,0.0) 300 continue lasl=lin lin=1 li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,3) do 40 ii=2,lasl lin=ii li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,2) if(mod(ii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 40 continue call linee c call newrgb(0.0,0.0,0.0) 100 CONTINUE call newrgb(0.0,0.0,0.0) 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) 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) call linee 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=XXL/FLOAT(NX1) DY=YYL/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 X1=0.5*DX*(2*I-NX2-1)+HX*FLOAT(NXP) Y1=0.5*DY*(2*J-NY2-1) Z1=0.0 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(M3) 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) C CALL PLOT(GX1,GY1,3) DO 20 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 R2=X2*X2+Y2*Y2+Z2*Z2 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 100 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 100 IF(R2.LT.AR2) GO TO 100 IF(Z2.GT.GZMA) GO TO 100 IF(Z2.LT.GZMI) GO TO 110 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 22 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) R2=X3*X3+Y3*Y3+Z3*Z3 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 100 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 100 IF(R2.LT.AR2) GO TO 100 IF(Z3.GT.GZMA) GO TO 100 IF(Z3.LT.GZMI) GO TO 110 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 22 CONTINUE C GX1=GX0+XDD1+AMP*(X3+AMPY*Y3*XC1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) C CALL PLOT(GX1,GY1,2) X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 24 II=1,8 Q(II)=Q(II+8) 24 CONTINUE 20 CONTINUE C GO TO 100 110 CONTINUE X1=0.5*DX*(2*I-NX2-1)+HX*FLOAT(NXP) Y1=0.5*DY*(2*J-NY2-1) Z1=0.0 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(M3) 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 if(mod(iii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 30 CONTINUE 100 CONTINUE RETURN END SUBROUTINE AINTEA(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 21 22 23 C GYMA DXL DYL 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) GYMA=AA(21) DXL=AA(22) DYL=AA(23) C XX=999.0 XM=0.15 YM=0.15 X1=GX0 Y1=GY0-2.0*YM c CALL CHARA1(X1,Y1,0,0,-1,0,MI,GXMI,GXMA) C C 1 X Y X=GX0 Y=GY0+0.0*GYL CALL PLOT(X,Y,3) X=GX0 Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=GX0+GXL CALL PLOT(X,Y,2) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) X=GX0 CALL PLOT(X,Y,2) C XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+0.5*GYL CALL PLOT(X,Y,3) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) C C 2 X Z X=GX0 Y=GY0+0.5*GYL+DYL CALL PLOT(X,Y,3) X=GX0 Y=GY0+1.0*GYL+DYL CALL PLOT(X,Y,2) X=GX0+GXL CALL PLOT(X,Y,2) Y=GY0+0.5*GYL+DYL CALL PLOT(X,Y,2) X=GX0 CALL PLOT(X,Y,2) C XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+1.0*GYL+DYL CALL PLOT(X,Y,3) Y=GY0+0.5*GYL+DYL CALL PLOT(X,Y,2) C C 3 Z Y X=GX0+GXL+DXL Y=GY0+0.0*GYL CALL PLOT(X,Y,3) X=GX0+GXL+DXL Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=GX0+GXL+0.5*GYL+DXL CALL PLOT(X,Y,2) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) X=GX0+GXL+DXL CALL PLOT(X,Y,2) call linee 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 C* 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 QUANT2(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 C 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 QUANT2(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 QUANT2(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*X1 GY1=GY0+0.5*GYL-AMP*Y1 CALL PLOT(GX1,GY1,3) C GX1=GX0+XDD1+AMP*X3 GY1=GY0+0.5*GYL-AMP*Y3 CALL PLOT(GX1,GY1,2) C GX2=GX0+XDD1+AMP*X1 GY2=GY0+0.0*GYL+AMP*Z1 GY2=0.5*GYL+DYL+GY2 CALL PLOT(GX2,GY2,3) C GX2=GX0+XDD1+AMP*X3 GY2=GY0+0.0*GYL+AMP*Z3 GY2=0.5*GYL+DYL+GY2 CALL PLOT(GX2,GY2,2) C GX3=GX0+AMP*Z1 GX3=GXL+DXL+GX3 GY3=GY0+0.5*GYL-AMP*Y1 CALL PLOT(GX3,GY3,3) C GX3=GX0+AMP*Z3 GX3=GXL+DXL+GX3 GY3=GY0+0.5*GYL-AMP*Y3 CALL PLOT(GX3,GY3,2) C X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE if(mod(iii,100).eq.0) call linee 30 CONTINUE 100 CONTINUE RETURN END SUBROUTINE AINTEB(IA,AA,F,P) PARAMETER(nline=2000,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO IC 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) IC=IA(9) 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) GYMA=AA(21) C XX=999.0 XM=0.15 YM=0.15 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 Y=GY0+0.0*GYL CALL PLOT(X,Y,3) if(ic.eq.2) go to 3 if(ic.eq.3) go to 4 if(ic.eq.4) go to 5 X=GX0 Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=GX0+GXL IF(IC.EQ.3) X=GX0+0.5*GYL CALL PLOT(X,Y,2) Y=GY0+0.0*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 C* CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1+GXL C* CALL PLOT(X,Y,2) X=GX0+GXL Y=GY0+0.5*GYL C* CALL PLOT(X,Y,2) C c IF(IC.EQ.3) GO TO 3 XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+0.5*GYL CALL PLOT(X,Y,3) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) X=X+0.5*GYL*XC1 Y=Y+0.5*GYL*XS1 C* CALL PLOT(X,Y,2) 3 CONTINUE if(ic.eq.1) go to 6 if(ic.eq.3) go to 4 if(ic.eq.4) go to 5 X=GX0 Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=GX0+GXL c IF(IC.EQ.3) X=GX0+0.5*GYL CALL PLOT(X,Y,2) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) X=GX0 CALL PLOT(X,Y,2) y=gy0-0.5*gyl call plot(x,y,2) x=gx0+gxl call plot(x,y,2) y=gy0 call plot(x,y,2) x=gx0 call plot(x,y,2) c11 X=GX0+0.5*GYL*XC1 c Y=GY0+0.5*GYL*XS1+0.5*GYL C* CALL PLOT(X,Y,2) c X=GX0+0.5*GYL*XC1+GXL C* CALL PLOT(X,Y,2) c X=GX0+GXL c Y=GY0+0.5*GYL C* CALL PLOT(X,Y,2) C c IF(IC.EQ.3) GO TO 3 XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+0.5*GYL CALL PLOT(X,Y,3) Y=GY0-0.5*GYL CALL PLOT(X,Y,2) c X=X+0.5*GYL*XC1 c Y=Y+0.5*GYL*XS1 C* CALL PLOT(X,Y,2) 4 continue if(ic.eq.1) go to 6 if(ic.eq.2) go to 6 if(ic.eq.4) go to 5 X=GX0 Y=GY0+0.0*GYL CALL PLOT(X,Y,3) c if(ic.eq.2) go to 3 c if(ic.eq.3) go to 4 X=GX0 Y=GY0+0.5*GYL CALL PLOT(X,Y,2) cc X=GX0+GXL X=GX0+0.5*GYL CALL PLOT(X,Y,2) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) X=GX0 CALL PLOT(X,Y,2) y=gy0+1.0*gyl call plot(x,y,2) x=gx0+0.5*gyl call plot(x,y,2) y=gy0+0.5*gyl call plot(x,y,2) x=gx0 call plot(x,y,2) 5 continue if(ic.eq.1) go to 6 if(ic.eq.2) go to 6 if(ic.eq.3) go to 6 c X=GX0 Y=GY0+0.0*GYL CALL PLOT(X,Y,3) c if(ic.eq.2) go to 3 c if(ic.eq.3) go to 4 X=GX0 Y=GY0+0.5*GYL CALL PLOT(X,Y,2) cc X=GX0+GXL X=GX0-0.5*GYL CALL PLOT(X,Y,2) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) X=GX0 CALL PLOT(X,Y,2) y=gy0+1.0*gyl call plot(x,y,2) x=gx0-0.5*gyl call plot(x,y,2) y=gy0+0.5*gyl call plot(x,y,2) x=gx0 call plot(x,y,2) 6 continue call linee 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 C* GYMA=0.5*(GXMA-GXMI) c GZMI=-0.1*HZ gymi=-gyma gzmi=-gyma 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 jj=1,2 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 if(jj.eq.1) Z1=SQRT(Z1) if(jj.eq.2) Z1=-SQRT(Z1) c Z1=SQRT(Z1) CALL QUANT2(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 IF(IC.NE.1) GO TO 11 GX1=GX0+XDD1+AMP*X1 GY1=GY0+0.5*GYL-AMP*Y1 11 CONTINUE IF(IC.NE.2) GO TO 12 GX1=GX0+XDD1+AMP*X1 GY1=GY0+0.0*GYL+AMP*Z1 12 CONTINUE IF(IC.NE.3) GO TO 13 c GX1=GX0+AMP*Z1 c GY1=GY0+0.5*GYL-AMP*Y1 GX1=GX0+AMP*y1 GY1=GY0+0.5*GYL+AMP*z1 13 CONTINUE IF(IC.NE.4) GO TO 14 c GX1=GX0+AMP*Z1 c GY1=GY0+0.5*GYL-AMP*Y1 GX1=GX0-AMP*y1 GY1=GY0+0.5*GYL-AMP*z1 14 CONTINUE CALL PLOT(GX1,GY1,3) lin=1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 c**** 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 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 CALL QUANT2(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) R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 CALL QUANT2(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 IF(IC.NE.1) GO TO 21 GX1=GX0+XDD1+AMP*X3 GY1=GY0+0.5*GYL-AMP*Y3 21 CONTINUE IF(IC.NE.2) GO TO 22 GX1=GX0+XDD1+AMP*X3 GY1=GY0+0.0*GYL+AMP*Z3 22 CONTINUE IF(IC.NE.3) GO TO 23 c GX1=GX0+AMP*Z3 c GY1=GY0+0.5*GYL-AMP*Y3 GX1=GX0+AMP*y3 GY1=GY0+0.5*GYL+AMP*z3 23 CONTINUE c IF(IC.NE.4) GO TO 24 c GX1=GX0+AMP*Z3 c GY1=GY0+0.5*GYL-AMP*Y3 GX1=GX0-AMP*y3 GY1=GY0+0.5*GYL-AMP*z3 24 CONTINUE c c CALL PLOT(GX1,GY1,2) lin=iii+1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,100).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c**** 200 continue if(lin.eq.1) go to 100 call newrgb(0.0,0.5,1.0) go to 300 210 continue if(lin.eq.1) go to 100 call newrgb(0.0,1.0,0.0) 300 continue lasl=lin lin=1 li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,3) do 40 ii=2,lasl lin=ii li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,2) if(mod(ii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 40 continue call linee 100 CONTINUE call newrgb(0.0,0.0,0.0) RETURN END SUBROUTINE AINTEC(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 IC 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) IC=IA(9) 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) GYMA=AA(21) C XX=999.0 XM=0.15 YM=0.15 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 Y=GY0+0.0*GYL CALL PLOT(X,Y,3) X=GX0 Y=GY0+0.5*GYL CALL PLOT(X,Y,2) X=GX0+GXL IF(IC.EQ.3) X=GX0+0.5*GYL CALL PLOT(X,Y,2) Y=GY0+0.0*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 C* CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1+GXL C* CALL PLOT(X,Y,2) X=GX0+GXL Y=GY0+0.5*GYL C* CALL PLOT(X,Y,2) C IF(IC.EQ.3) GO TO 3 XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+0.5*GYL CALL PLOT(X,Y,3) Y=GY0+0.0*GYL CALL PLOT(X,Y,2) X=X+0.5*GYL*XC1 Y=Y+0.5*GYL*XS1 C* CALL PLOT(X,Y,2) 3 CONTINUE call linee 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 C* 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 C JB1=NY2/2+1 C DO 100 J=JB1,NY2 C DO 100 I=1,NX2 C X=0.5*DX*(2*I-NX2-1) C Y=0.5*DY*(2*J-NY2-1) DO 100 J=1,NY1 DO 100 I=1,NX TH1=PI*FLOAT(J-1)/FLOAT(NY) TH2=PI*(0.5-FLOAT(I)/(6.0*FLOAT(NX))) X=COS(TH2)*COS(TH1) Y=COS(TH2)*SIN(TH1) C 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 QUANT2(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 IF(IC.NE.1) GO TO 11 GX1=GX0+XDD1+AMP*X1 GY1=GY0+0.5*GYL-AMP*Y1 11 CONTINUE IF(IC.NE.2) GO TO 12 GX1=GX0+XDD1+AMP*X1 GY1=GY0+0.0*GYL+AMP*Z1 12 CONTINUE IF(IC.NE.3) GO TO 13 GX1=GX0+AMP*Z1 GY1=GY0+0.5*GYL-AMP*Y1 13 CONTINUE 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 QUANT2(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.0.0) GO TO 40 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 100 CALL QUANT2(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 GO TO 42 40 CONTINUE Z2=Z3 Z3=0.0 X3=X1+(X3-X1)*(Z3-Z1)/(Z2-Z1) Y3=Y1+(Y3-Y1)*(Z3-Z1)/(Z2-Z1) 42 CONTINUE C IF(IC.NE.1) GO TO 21 GX1=GX0+XDD1+AMP*X3 GY1=GY0+0.5*GYL-AMP*Y3 21 CONTINUE IF(IC.NE.2) GO TO 22 GX1=GX0+XDD1+AMP*X3 GY1=GY0+0.0*GYL+AMP*Z3 22 CONTINUE IF(IC.NE.3) GO TO 23 GX1=GX0+AMP*Z3 GY1=GY0+0.5*GYL-AMP*Y3 23 CONTINUE CALL PLOT(GX1,GY1,2) IF(Z2.LT.0.0) GO TO 100 X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE if(mod(iii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 30 CONTINUE 100 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 FILE NAME / POSC / SUBPS3 K. OHTA 03/17/1992 c file name /test3/subc.f c PostScript Program File. c Created by T. Ogino and K. Ohta at April 1992 c STELab Nagoya University. subroutine xyopen write(10,10) 10 format( "%!" ) write(10,20) 20 format( "%%Page: 1 1" ) write(10,*) "/m { moveto } def" write(10,*) "/l { lineto } def" write(10,*) "/s { show } def" write(10,*) "/aws { awidthshow } def" write(10,*) "/sl { setlinewidth } def" write(10,*) "/sn { stroke newpath } def" write(10,*) "/tr { /Times-Roman findfont } def" write(10,*) "/sb { /Symbol findfont } def" write(10,*) "/sf { scalefont } def" write(10,*) "/se { setfont } def" write(10,*) "/ro { rotate } def" write(10,*) "/sd { setdash } def" write(10,*) "/st { stroke } def" write(10,*) "/sp { stroke showpage } def" write(10,*) "/tl { translate } def" write(10,*) "/sc { scale } def" write(10,*) "/np { newpath } def" write(10,*) "/cp { closepath } def" write(10,*) "/gs { gsave } def" write(10,*) "/gr { grestore } def" write(10,*) "/grc { grestore 0.0 0.0 0.0 setrgbcolor } def" write(10,*) "/sg { setgray } def" write(10,*) "/chs { sethsbcolor } def" write(10,*) "/crg { setrgbcolor } def" write(10,*) "/f { fill } def" write(10,*) "/cup { currentpoint } def" write(10,*) "/ctr { /Times-Roman findfont } def" write(10,*) "/cti { /Times-Italic findfont } def" write(10,*) "/ctb { /Times-Bold findfont } def" write(10,*) "/ctbi { /Times-BoldItalic findfont } def" write(10,*) "/chr { /Helvetica findfont } def" write(10,*) "/cho { /Helvetica-Oblique findfont } def" write(10,*) "/chb { /Helvetica-Bold findfont } def" write(10,*) "/chbo { /Helvetica-BoldOblique findfont } def" write(10,*) "/ccr { /Courier findfont } def" write(10,*) "/cco { /Courier-Oblique findfont } def" write(10,*) "/ccb { /Courier-Bold findfont } def" write(10,*) "/ccbo { /Courier-BoldOblique findfont } def" write(10,*) "/c {setrgbcolor 1 0 translate newpath 0 0 moveto & 0 1 lineto 1 1 lineto 1 0 lineto closepath fill} def" write(10,*) "/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(10,10) ino,ino 10 format("%%Page:",1x,i2,1x,i2) c if(idci.eq.1) then if(ity.eq.1.and.(isz.eq.1.or.isz.eq.2)) then write(10,*) " 1.00 1.00 sc" write(10,*) " 28.0 33.0 tl" else if(ity.eq.2.and.isz.eq.1) then write(10,*) " 1.00 1.00 sc" write(10,*) " 90.0 ro" write(10,*) " 35.0 -565.0 tl" else if(ity.eq.2.and.isz.eq.2) then write(10,*) " 1.00 1.00 sc" write(10,*) " 90.0 ro" write(10,*) " 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(10,*) " 28.38 28.17 sc" write(10,*) " 1.00 1.20 tl" else if(ity.eq.2.and.isz.eq.1) then write(10,*) " 90.0 ro" write(10,*) " 28.17 28.38 sc" write(10,*) " 1.15 -19.90 tl" else if(ity.eq.2.and.isz.eq.2) then write(10,*) " 90.0 ro" write(10,*) " 28.17 28.38 sc" write(10,*) " 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(10,*) " 72.08 71.55 sc" write(10,*) " 0.39 0.47 tl" else if(ity.eq.2.and.isz.eq.1) then write(10,*) " 90.0 ro" write(10,*) " 71.55 72.08 sc" write(10,*) " 0.49 -7.88 tl" else if(ity.eq.2.and.isz.eq.2) then write(10,*) " 90.0 ro" write(10,*) " 71.55 72.08 sc" write(10,*) " 0.49 -9.66 tl" end if end if call newpen(1) return end subroutine plote write(10,*) "sp" return end subroutine linee write(10,*) "st" return end subroutine dashe write(10,*) "sn [ ] 0 sd" return end subroutine plot(x1,y1,im) common /cdci/ldci if(ldci.eq.1) then if(im.eq.3) then write(10,10) x1,y1 10 format(f6.1,1x,f6.1," m") else if(im.eq.2) then write(10,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(10,30) x1,y1 30 format(f6.2,1x,f6.2," m") else if(im.eq.2) then write(10,40) x1,y1 40 format(f6.2,1x,f6.2," l") end if end if return end subroutine factor(fct) write(10,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(10,30) fp 30 format("sn",1x,f6.3," sl") if(ipn.gt.10) then if(ldci.eq.1) then if(is.eq.1) write(10,*) "[50.0 0.00] 0 sd" if(is.eq.2) write(10,*) "[1.50 1.50] 0 sd" if(is.eq.3) write(10,*) "[3.00 3.00] 0 sd" if(is.eq.4) write(10,*) "[6.00 6.00] 0 sd" if(is.eq.5) write(10,*) "[9.00 9.00] 0 sd" if(is.eq.6) write(10,*) "[12.0 12.0] 0 sd" if(is.eq.7) write(10,*) "[6.00 3.00] 0 sd" if(is.eq.8) write(10,*) "[12.0 3.00] 0 sd" if(is.eq.9) write(10,*) "[18.0 3.00] 0 sd" if(is.eq.10) write(10,*) "[24.0 3.00] 0 sd" else if(ldci.eq.2) then if(is.eq.1) write(10,*) "[1.00 0.00] 0 sd" if(is.eq.2) write(10,*) "[0.05 0.05] 0 sd" if(is.eq.3) write(10,*) "[0.10 0.10] 0 sd" if(is.eq.4) write(10,*) "[0.20 0.20] 0 sd" if(is.eq.5) write(10,*) "[0.30 0.30] 0 sd" if(is.eq.6) write(10,*) "[0.40 0.40] 0 sd" if(is.eq.7) write(10,*) "[0.20 0.10] 0 sd" if(is.eq.8) write(10,*) "[0.40 0.10] 0 sd" if(is.eq.9) write(10,*) "[0.60 0.10] 0 sd" if(is.eq.10) write(10,*) "[0.80 0.10] 0 sd" else if(ldci.eq.3) then if(is.eq.1) write(10,*) "[1.00 0.00] 0 sd" if(is.eq.2) write(10,*) "[0.02 0.02] 0 sd" if(is.eq.3) write(10,*) "[0.04 0.04] 0 sd" if(is.eq.4) write(10,*) "[0.08 0.08] 0 sd" if(is.eq.5) write(10,*) "[0.12 0.12] 0 sd" if(is.eq.6) write(10,*) "[0.16 0.16] 0 sd" if(is.eq.7) write(10,*) "[0.08 0.04] 0 sd" if(is.eq.8) write(10,*) "[0.16 0.04] 0 sd" if(is.eq.9) write(10,*) "[0.24 0.04] 0 sd" if(is.eq.10) write(10,*) "[0.32 0.04] 0 sd" end if end if return end subroutine symbol(x,y,h,isymb,ang,n) character isymb*80,ica*80,ich(80)*1 equivalence (ica,ich) ica=isymb write(10,*) "gs" h0=h*1.50 write(10,10) h0 10 format("tr ",f6.2," sf se") if(x.gt.998.0.or.y.gt.998.0) then write(10,*) "cup m" else write(10,20) x,y 20 format(f6.2,1x,f6.2," m") end if write(10,30) ang 30 format(f6.1," ro") h1=h*0.08 h2=h*0.06 write(10,40) h1,h2 40 format(f6.3," 0.0 8#040 ",f6.3," 0.0") c write(10,*) "(",(ich(i),i=1,n),") s" write(10,*) "(",(ich(i),i=1,n),") aws" write(10,*) "gr" return end subroutine ctimes(ichr) common /cchr/lchr lchr=min0(ichr,4) lchr=max0(ichr,1) return end subroutine chelve(ichr) common /cchr/lchr lchr=min0(ichr,4)+4 lchr=max0(ichr,1)+4 return end subroutine ccouri(ichr) common /cchr/lchr lchr=min0(ichr,4)+8 lchr=max0(ichr,1)+8 return end subroutine symblb(x,y,h,isymb,ang,n) character isymb*80,ica*80,ich(80)*1,icb(12)*4 equivalence (ica,ich) common /cchr/lchr data icb/"ctr ","cti ","ctb ","ctbi","chr ","cho ", & "chb ","chbo","ccr ","cco ","ccb ","ccbo"/ ica=isymb write(10,*) "gs" if(lchr.eq.0) write(10,*) "tr" if(lchr.ge.1) write(10,*) icb(lchr) h0=h*1.50 write(10,10) h0 10 format(f5.2," sf se") if(x.gt.998.0.or.y.gt.998.0) then write(10,*) "cup m" else write(10,20) x,y 20 format(f6.2,1x,f6.2," m") end if write(10,30) ang 30 format(f6.1," ro") if(lchr.le.8) then h1=h*0.08 h2=h*0.06 write(10,40) h1,h2 40 format(f6.3," 0.0 8#040 ",f6.3," 0.0") write(10,*) "(",(ich(i),i=1,n),") aws" else if(lchr.ge.9) then write(10,*) "(",(ich(i),i=1,n),") s" end if write(10,*) "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(10,*) "gs" h0=h*1.50 write(10,10) h0 10 format("sb ",f6.2," sf se") if(x.gt.998.0.or.y.gt.998.0) then write(10,*) "cup m" else write(10,20) x,y 20 format(f6.2,1x,f6.2," m") end if write(10,30) ang 30 format(f6.1," ro") h1=h*0.08 h2=h*0.06 write(10,40) h1,h2 40 format(f6.3," 0.0 8#040 ",f6.3," 0.0") c write(10,*) "(",(ich(i),i=1,n),") s" write(10,*) "(",(ich(i),i=1,n),") aws" write(10,*) "gr" return end subroutine circ1(x,y,r) x1=x+r write(10,10) x1,y 10 format(f6.2,1x,f6.2," m") write(10,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(10,*) "gs" return end subroutine grae write(10,*) "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(10,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(10,10) ar,ag,ab 10 format(3(1x,f4.2)," crg") return end subroutine newgry(ig) is=max0(ig,0) is=min0(ig,15) fg=float(is)*1.0/15.0 write(10,10) fg 10 format(f5.3," sg") return end subroutine paint write(10,*) "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(10,*) "gs" write(10,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(10,22) xa,ya c do 100 j=1,ny y=1.0 x=-float(nx) write(10,22) x,y 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(10,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(10,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(10,*) "st" write(10,*) "grc" return end SUBROUTINE AINTE11(IA,AA,F,P) PARAMETER(nline=2000,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt) 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) COTL=AA(25) SITL=AA(26) c ccc waku kaku tokoro 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 y=gy0 call plot(x,y,2) X=GX0 CALL PLOT(X,Y,2) 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 X=GX0+xdd1 Y=GY0+GYL*0.5 CALL PLOT(X,Y,2) Y=GY0 CALL PLOT(X,Y,2) c call linee ccc waku kakiowatta! C ccc boundary?? condition ?? 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) c GZMI=-0.1*HZ gymi=-gyma gzmi=-gyma c GZMI GZMA ?? GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 C cc polar cap kara start suru basyo wo kimeru JB1=NY2/2+1 do 100 jj=1,2 DO 100 J=JB1,NY2 cc DO 100 J=1,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 c z ga + sika naiyouna.... Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=ARB*ARB-X1*X1-Y1*Y1 c zga 0.0 yori tiisaitoki 100 he IF(Z1.LE.0.0) GO TO 100 c IF(Z1.LE.0.0) GO TO 1030 if(jj.eq.1) Z1=SQRT(Z1) if(jj.eq.2) Z1=-SQRT(Z1) c X2=X1*COTL-Z1*SITL Z2=X1*SITL+Z1*COTL c IF(AA(25).LT.0.0) X2=-X2 c IF(AA(25).LT.0.0) Z2=-Z2 X1=X2 Z1=Z2 c 1040 CALL QUANT11(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) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c CALL PLOT(GX1,GY1,3) lin=1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 c**** 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 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 CALL QUANT11(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) R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 CALL QUANT11(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) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) c CALL PLOT(GX1,GY1,2) lin=iii+1 li1=2*(lin-1)+1 li2=2*lin po(li1)=GX1 po(li2)=GY1 X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,100).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c**** 200 continue if(lin.eq.1) go to 100 call newrgb(0.0,0.5,1.0) go to 300 210 continue if(lin.eq.1) go to 100 call newrgb(0.0,1.0,0.0) 300 continue lasl=lin lin=1 li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,3) do 40 ii=2,lasl lin=ii li1=2*(lin-1)+1 li2=2*lin GX1=po(li1) GY1=po(li2) call plot(gx1,gy1,2) if(mod(ii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 40 continue call linee 100 CONTINUE call newrgb(0.0,0.0,0.0) RETURN END SUBROUTINE QUANT11(X,Y,Z,HX,HY,HZ,MM,NX,NY,NZ,NXP,AA,F,P) DIMENSION F(1),P(1),AA(1),B(3) C xo=x yo=y zo=z if(yo.lt.0.0) then y=-y z=z endif NX1=NX+1 NX2=NX+2 NY1=NY+1 NZ1=NZ+1 NY2=NY+2 NZ2=NZ+2 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) COTL=AA(25) SITL=AA(26) C XI=X/HX+0.5*FLOAT(NX2+1-2*NXP) XJ=Y/HY+1.5 c XK=Z/HZ+0.5*FLOAT(NZ+3) 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,3 I1=I1+N3 I2=M+L1 P(I2)=XI*XJ*XK*F(I1+N22)+XI*XJ*XK1*F(I1+N11) 1 +XI*XJ1*XK*F(I1+N12)+XI1*XJ*XK*F(I1+N21) 2 +XI*XJ1*XK1*F(I1+1)+XI1*XJ*XK1*F(I1+N1) 3 +XI1*XJ1*XK*F(I1+N2)+XI1*XJ1*XK1*F(I1) P(I2+5)=P(I2) 10 CONTINUE C XT=X*COTL+Z*SITL ZT=-X*SITL+Z*COTL RO2=X*X+Y*Y+Z*Z RO1=SQRT(RO2) RO3=RO1*RO2 RO5=RO2*RO3 AR2=AR1*AR1 X2=0.0 IF(RO2.GT.AR2) X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) BX=-3.0*B0*XT*ZT/RO5 BY=-3.0*B0*Y*ZT/RO5 BZ=B0*(XT*XT+Y*Y-2.0*ZT*ZT)/RO5 B(1)=BX*COTL-BZ*SITL B(2)=BY B(3)=BX*SITL+BZ*COTL 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 if(yo.lt.0.0) then i1=2+L1 p(i1)=-p(i1) endif x=xo y=yo z=zo RETURN END SUBROUTINE AINTE21(IA,AA,F,P) PARAMETER(nline=2000,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt,2),ipo(3,2) 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) IPRO=IA(10) 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) XLIM=AA(23) XLI2=AA(24) 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 C y=gy0 call plot(x,y,2) X=GX0 CALL PLOT(X,Y,2) 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 X=GX0+xdd1 Y=GY0+GYL*0.5 CALL PLOT(X,Y,2) Y=GY0 CALL PLOT(X,Y,2) c call linee 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=XXL/FLOAT(NX1) c DY=YYL/FLOAT(NY1) DY=2.0*YYL/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) gymi=-gyma AMPY=0.5*(GXMA-GXMI)*GYL/(GYMA*GXL) GZMI=-0.1*HZ c GZMI=-GYMA 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 c do 100 jj=1,2 DO 100 J=JB1,NY2 C DO 100 J=1,NY1 c DO 100 J=1,NY2 DO 100 I=1,NX2 c do 300 jj=1,2 do 300 jj=1,1 X1=0.5*DX*(2*I-NX2-1)+HX*FLOAT(NXP) Y1=0.5*DY*(2*J-NY2-1) C Y1=0.5*DY*(J) C y1=0.5*(2*J-1) Z1=0.0 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 300 AA(11)=1.0 X2=Q(M3) IF(X2.LT.0.0) AA(11)=-1.0 if(jj.eq.2) aa(11)=-aa(11) 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) IF(IPRO.EQ.1) GX1=GX0+XDD1+AMP*X1 IF(IPRO.EQ.1) GY1=GY0+0.5*GYL+AMP*AMPY*Z1 IF(IPRO.EQ.2) GX1=GX0+XDD1+AMP*X1 IF(IPRO.EQ.2) GY1=GY0+0.5*GYL-AMP*AMPY*Y1 c CALL PLOT(GX1,GY1,3) lin=1 li1=2*(lin-1)+1 li2=2*lin po(li1,jj)=GX1 po(li2,jj)=GY1 c**** 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 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 c XLIM=0.0 IF(X2.LT.XLIM.AND.Y2.LT.GYMI) GO TO 300 IF(X2.LT.XLIM.AND.Y2.GT.GYMA) GO TO 300 IF(X2.LT.XLIM.AND.Z2.LT.GZMI) GO TO 300 IF(X2.LT.XLIM.AND.Z2.GT.GZMA) GO TO 300 IF(X2.LT.XLI2.AND.Y2.LT.GYMI) GO TO 220 IF(X2.LT.XLI2.AND.Y2.GT.GYMA) GO TO 220 IF(X2.LT.XLI2.AND.Z2.LT.GZMI) GO TO 240 IF(X2.LT.XLI2.AND.Z2.GT.GZMA) GO TO 220 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 c IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 IF(Z2.GT.GZMA) GO TO 200 IF(Z2.LT.GZMI) GO TO 240 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 300 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) R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.XLIM.AND.Y3.LT.GYMI) GO TO 300 IF(X3.LT.XLIM.AND.Y3.GT.GYMA) GO TO 300 IF(X3.LT.XLIM.AND.Z3.LT.GZMI) GO TO 300 IF(X3.LT.XLIM.AND.Z3.GT.GZMA) GO TO 300 IF(X3.LT.XLI2.AND.Y3.LT.GYMI) GO TO 220 IF(X3.LT.XLI2.AND.Y3.GT.GYMA) GO TO 220 IF(X3.LT.XLI2.AND.Z3.LT.GZMI) GO TO 240 IF(X3.LT.XLI2.AND.Z3.GT.GZMA) GO TO 220 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 c IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 IF(Z3.GT.GZMA) GO TO 200 IF(Z3.LT.GZMI) GO TO 240 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 300 32 CONTINUE C GX1=GX0+XDD1+AMP*(X3+AMPY*Y3*XC1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) IF(IPRO.EQ.1) GX1=GX0+XDD1+AMP*X3 IF(IPRO.EQ.1) GY1=GY0+0.5*GYL+AMP*AMPY*Z3 IF(IPRO.EQ.2) GX1=GX0+XDD1+AMP*X3 IF(IPRO.EQ.2) GY1=GY0+0.5*GYL-AMP*AMPY*Y3 c CALL PLOT(GX1,GY1,2) lin=iii+1 li1=2*(lin-1)+1 li2=2*lin po(li1,jj)=GX1 po(li2,jj)=GY1 X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,300).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c go to 300 c**** 200 continue ipo(1,jj)=lin ipo(2,jj)=4 c if(lin.eq.1) go to 100 c call newrgb(1.0,0.0,0.0) go to 300 220 continue ipo(1,jj)=lin ipo(2,jj)=2 c if(lin.eq.1) go to 100 c call newrgb(1.0,0.8,0.0) go to 300 210 continue ipo(1,jj)=lin ipo(2,jj)=3 c if(lin.eq.1) go to 100 c call newrgb(0.0,1.0,0.0) go to 300 240 continue ipo(1,jj)=lin ipo(2,jj)=1 c if(lin.eq.1) go to 100 c call newrgb(1.0,0.0,0.0) 300 continue c line1=ipo(1,1) line2=ipo(1,2) col1=ipo(2,1) col2=ipo(2,2) line1=min(line1,line2) c if(line1.eq.1) go to 100 c if(col1.eq.1.and.col2.eq.1) color=1 c if(col1.eq.1.and.col2.eq.2) color=2 c if(col1.eq.1.and.col2.eq.3) color=4 c if(col1.eq.2.and.col2.eq.1) color=2 c if(col1.eq.2.and.col2.eq.2) color=2 c if(col1.eq.2.and.col2.eq.3) color=4 c if(col1.eq.3.and.col2.eq.1) color=4 c if(col1.eq.3.and.col2.eq.2) color=4 c if(col1.eq.3.and.col2.eq.3) color=3 c ipo(2,1)=color c ipo(2,2)=color c do 42 jj=1,1 lin=ipo(1,jj) if(lin.eq.1) go to 42 icol=ipo(2,jj) if(icol.eq.1) call newrgb(1.0,0.0,0.0) if(icol.eq.2) call newrgb(1.0,0.8,0.0) if(icol.eq.3) call newrgb(0.0,1.0,0.0) if(icol.eq.4) call newrgb(1.0,0.0,1.0) c lasl=lin lin=1 li1=2*(lin-1)+1 li2=2*lin GX1=po(li1,jj) GY1=po(li2,jj) call plot(gx1,gy1,3) do 40 ii=2,lasl lin=ii li1=2*(lin-1)+1 li2=2*lin GX1=po(li1,jj) GY1=po(li2,jj) call plot(gx1,gy1,2) if(mod(ii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 40 continue call linee 42 continue c call newrgb(0.0,0.0,0.0) 100 CONTINUE call newrgb(0.0,0.0,0.0) RETURN END SUBROUTINE AINTE22(IA,AA,F,P) PARAMETER(nline=2000,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt,2),ipo(3,2) 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) XLIM=AA(23) XLI2=AA(24) 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 C y=gy0 call plot(x,y,2) X=GX0 CALL PLOT(X,Y,2) 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 X=GX0+xdd1 Y=GY0+GYL*0.5 CALL PLOT(X,Y,2) Y=GY0 CALL PLOT(X,Y,2) c call linee 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=XXL/FLOAT(NX1) DY=YYL/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) gymi=-gyma AMPY=0.5*(GXMA-GXMI)*GYL/(GYMA*GXL) c GZMI=-0.1*HZ GZMI=-GYMA 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 c do 100 jj=1,2 DO 100 J=JB1,NY2 C DO 100 J=1,NY1 c DO 100 J=1,NY2 DO 100 I=1,NX2 do 300 jj=1,2 X1=0.5*DX*(2*I-NX2-1)+HX*FLOAT(NXP) Y1=0.5*DY*(2*J-NY2-1) C Y1=0.5*DY*(J) C y1=0.5*(2*J-1) Z1=0.0 CALL QUANT11(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 300 AA(11)=1.0 X2=Q(M3) IF(X2.LT.0.0) AA(11)=-1.0 if(jj.eq.2) aa(11)=-aa(11) 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) c CALL PLOT(GX1,GY1,3) lin=1 li1=2*(lin-1)+1 li2=2*lin po(li1,jj)=GX1 po(li2,jj)=GY1 c**** 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 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 c XLIM=0.0 IF(X2.LT.XLIM.AND.Y2.LT.GYMI) GO TO 300 IF(X2.LT.XLIM.AND.Y2.GT.GYMA) GO TO 300 IF(X2.LT.XLIM.AND.Z2.LT.GZMI) GO TO 300 IF(X2.LT.XLIM.AND.Z2.GT.GZMA) GO TO 300 IF(X2.LT.XLI2.AND.Y2.LT.GYMI) GO TO 220 IF(X2.LT.XLI2.AND.Y2.GT.GYMA) GO TO 220 IF(X2.LT.XLI2.AND.Z2.LT.GZMI) GO TO 220 IF(X2.LT.XLI2.AND.Z2.GT.GZMA) GO TO 220 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 CALL QUANT11(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 300 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) R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.XLIM.AND.Y3.LT.GYMI) GO TO 300 IF(X3.LT.XLIM.AND.Y3.GT.GYMA) GO TO 300 IF(X3.LT.XLIM.AND.Z3.LT.GZMI) GO TO 300 IF(X3.LT.XLIM.AND.Z3.GT.GZMA) GO TO 300 IF(X3.LT.XLI2.AND.Y3.LT.GYMI) GO TO 220 IF(X3.LT.XLI2.AND.Y3.GT.GYMA) GO TO 220 IF(X3.LT.XLI2.AND.Z3.LT.GZMI) GO TO 220 IF(X3.LT.XLI2.AND.Z3.GT.GZMA) GO TO 220 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 CALL QUANT11(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 300 32 CONTINUE C GX1=GX0+XDD1+AMP*(X3+AMPY*Y3*XC1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) c CALL PLOT(GX1,GY1,2) lin=iii+1 li1=2*(lin-1)+1 li2=2*lin po(li1,jj)=GX1 po(li2,jj)=GY1 X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,300).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c go to 300 c**** 200 continue ipo(1,jj)=lin ipo(2,jj)=1 c if(lin.eq.1) go to 100 c call newrgb(1.0,0.0,0.0) go to 300 220 continue ipo(1,jj)=lin ipo(2,jj)=2 c if(lin.eq.1) go to 100 c call newrgb(1.0,0.8,0.0) go to 300 210 continue ipo(1,jj)=lin ipo(2,jj)=3 c if(lin.eq.1) go to 100 c call newrgb(0.0,1.0,0.0) 300 continue c line1=ipo(1,1) line2=ipo(1,2) col1=ipo(2,1) col2=ipo(2,2) line1=min(line1,line2) if(line1.eq.1) go to 100 if(col1.eq.1.and.col2.eq.1) color=1 if(col1.eq.1.and.col2.eq.2) color=2 if(col1.eq.1.and.col2.eq.3) color=4 if(col1.eq.2.and.col2.eq.1) color=2 if(col1.eq.2.and.col2.eq.2) color=2 if(col1.eq.2.and.col2.eq.3) color=4 if(col1.eq.3.and.col2.eq.1) color=4 if(col1.eq.3.and.col2.eq.2) color=4 if(col1.eq.3.and.col2.eq.3) color=3 ipo(2,1)=color ipo(2,2)=color c do 42 jj=1,2 lin=ipo(1,jj) if(lin.eq.1) go to 42 icol=ipo(2,jj) if(icol.eq.1) call newrgb(1.0,0.0,0.0) if(icol.eq.2) call newrgb(1.0,0.8,0.0) if(icol.eq.3) call newrgb(0.0,1.0,0.0) if(icol.eq.4) call newrgb(0.0,0.5,1.0) c lasl=lin lin=1 li1=2*(lin-1)+1 li2=2*lin GX1=po(li1,jj) GY1=po(li2,jj) call plot(gx1,gy1,3) do 40 ii=2,lasl lin=ii li1=2*(lin-1)+1 li2=2*lin GX1=po(li1,jj) GY1=po(li2,jj) call plot(gx1,gy1,2) if(mod(ii,100).eq.0) then call linee call plot(gx1,gy1,3) end if 40 continue call linee 42 continue c call newrgb(0.0,0.0,0.0) 100 CONTINUE call newrgb(0.0,0.0,0.0) RETURN END