C SUBROUTINE TPERSP ====*====3====*====4====*====5====*====6====*====7
C
C THREE DIMENSIONAL PLOT PROGRAM SET
C
C 1989/12/09 BY A.YANASE
C 2001/03/09 BY A.YANASE
C THIS FILE CONTAINS THE FUNCTIONS
C ONLY FOR TPFXYZ & TPSECT
C 2004/05/10 MODIFIED BY H.FUNASHIMA
C BUG FIX FOR
C SYNTAX ERROR ABOUT DOUBLE PRECISION
C IN SOME SUBROUTINE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPERSP(XFU,XFG,ALP,BET,GAM,ELD,ELL)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW,IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
COMMON/TPERS6/DFST,DMIN,DMAX,DFAC,DVID
COMMON/TPERS6D/DDFST,DDMIN,DDMAX,DDFAC,DDVID
COMMON/TPERS9/XOM(3),WM(3,2)
C INTEGER TPNCUT
C EXTERNAL FUNC,TPDUMY,TPLTRA
C EXTERNAL TPNCUT
DIMENSION XFU(3,2),XFG(3)
WRITE(6,*) ' ----- WELCOME TO TPERSP V3.0 2001/07/07 -----'
IF(ALP.GE.180.0.OR.ALP.LE.-180.) GO TO 999
IF(BET.GE.90.0.OR.BET.LE.-90.0) GO TO 999
IF(GAM.GE.90.0.OR.GAM.LE.-90.0) GO TO 999
IF(ELD.LT.0.0) GO TO 999
IF(ELL.LT.0.0) GO TO 999
DO 10 I=1,3
IF(XFU(I,1).GE.XFU(I,2)) GO TO 999
IF(XFG(I).LE.0.0) GO TO 999
10 CONTINUE
C
IUOD=1
IHL=0
IPR=0
IMOD=1
CALL TPHILP(-1.0D0)
CALL TPHILD(0.0D0,0,0,0,0)
C ****
ALPHA=ALP
BETA=BET
GAMM=GAM
X0=XFU(1,1)
Y0=XFU(2,1)
Z0=XFU(3,1)
XN=XFU(1,2)
YN=XFU(2,2)
ZN=XFU(3,2)
XW=XFG(1)
YW=XFG(2)
ZW=XFG(3)
DO 1 I=1,3
XOM(I)=XFU(I,1)
DO 2 J=1,2
2 WM(I,J)=0.0
1 CONTINUE
WM(1,1)=1.0
WM(2,2)=1.0
XFCTR=XW/(XN-X0)
YFCTR=YW/(YN-Y0)
ZFCTR=ZW/(ZN-Z0)
C ****
RAD=3.141593/180.0
C ****
CA=COS(RAD*ALPHA)
SA=SIN(RAD*ALPHA)
XWC=XW*CA
XWS=XW*SA
YWC=YW*CA
YWS=YW*SA
C ****
IF(ALPHA.LT.-90.0) GO TO 301
IF(ALPHA.LT. 0.0) GO TO 302
IF(ALPHA.LT. 90.0) GO TO 303
GO TO 304
301 TAU=-(XWS+YWC)
GO TO 305
302 TAU=-XWS
GO TO 305
303 TAU=0.0
GO TO 305
304 TAU=-YWC
305 CONTINUE
EL=ELL
TAU=TAU+ELD
C ****
ET=EL+TAU
SG=ET*SIN(RAD*BETA)/COS(RAD*BETA)
C ****
V1=SG/ET
V2=(SG-YWS)/(ET+YWC)
V3=(SG+XWC)/(ET+XWS)
V4=(SG+XWC-YWS)/(ET+XWS+YWC)
XP0=MIN(V1,V2,V3,V4)*EL
C ****
VMIN=EL+ELD
TG=SIN(RAD*GAMM)/COS(RAD*GAMM)
ZC=VMIN*TG
ZM=TG
IF(GAMM.GT.0.0) ZM=VMIN/
& (VMIN+ABS(XWS)+ABS(YWC))*TG
C ****
XEYE=-(ET*SA+SG*CA)
YEYE=-(ET*CA-SG*SA)
ZEYE=-VMIN*TG
DFST=0.2
DMIN=0.1
DMAX=10.0
DFAC=2.0
DVID=0.25
DDFST=0.1
DDMIN=0.1
DDMAX=10.0
DDFAC=2.0
DDVID=0.25
WRITE(6,600) X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW
WRITE(6,603) XFCTR,YFCTR,ZFCTR
WRITE(6,601) ALPHA,BETA,GAMM,ELD,ELL
WRITE(6,602) XEYE,YEYE,ZEYE
600 FORMAT(' YOUR FIGURE IS FROM (',3E14.5,')'
& /' TO (',3E14.5,')'
& /' SIZE OF FIGURE (',3E14.5,')')
603 FORMAT(' SETTING FACTOR (',3E14.5,')')
601 FORMAT(' ALPH=',E14.5,' BETA=',E14.5,' GAMM=',E14.5,
& /' ELD= ',E14.5,' ELL= ',E14.5)
602 FORMAT(' EYE POINT IS (',3E14.5,')')
CALL TPPLOT(X0,Y0,Z0,3)
C this call is 1994/10/14
RETURN
999 CONTINUE
WRITE(6,690)
690 FORMAT(' *** INVALID VALUE OF ARGUMENT IN TPERSP ***')
RETURN
END
C SUBROUTINE TPSCHG ====*====3====*====4====*====5====*====6====*====7
C
C MODE OR UNIT CHANGE FOR TPERSP
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPSCHG(IA)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(11)
COMMON/TPERS6/DFST,DMIN,DMAX,DFAC,DVID
COMMON/TPERS9/XOM(3),WM(3,2)
DIMENSION AX(3)
IF(IABS(IA).EQ.1) IUOD=IA
RETURN
ENTRY TPCIHL(IA)
IF(IA.GE.-2.AND.IA.LE.5) IHL=IA
RETURN
ENTRY TPTRAC(IA)
IF(IA.LE.2.AND.IA.GE.0) IPR=IA
RETURN
ENTRY TPUNIT(IA)
IF(IA.EQ.IMOD) GO TO 11
IF(IABS(IA).NE.1) GO TO 11
IF(IA.EQ.1) GO TO 2
DO 3 J=1,2
DO 4 I=1,3
4 AX(I)=WM(I,J)*FCTR(I)
WA=SQRT(AX(1)*AX(1)+AX(2)*AX(2)+AX(3)*AX(3))
IF(WA.LT.1.0D-30) GO TO 1
DO 5 I=1,3
5 WM(I,J)=AX(I)/WA
3 CONTINUE
DO 6 I=1,3
6 XOM(I)=(XOM(I)-XM(I,1))*FCTR(I)
GO TO 1
2 DO 7 J=1,2
DO 8 I=1,3
8 AX(I)=WM(I,J)/FCTR(I)
WA=SQRT(AX(1)*AX(1)+AX(2)*AX(2)+AX(3)*AX(3))
IF(WA.LT.1.0D-30) GO TO 1
DO 9 I=1,3
9 WM(I,J)=AX(I)/WA
7 CONTINUE
DO 10 I=1,3
10 XOM(I)=XOM(I)/FCTR(I)+XM(I,1)
1 IMOD=IA
11 RETURN
ENTRY TPHLDS(DF,DI,DA,DC,DD)
IF(DF.GT.1.0D-10) DFST=DF
IF(DI.GT.1.0D-10) DMIN=DI
IF(DA.GT.1.0D-10) DMAX=DA
IF(DC.GT.1.0D-10) DFAC=DC
IF(DD.GT.1.0D-10) DVID=DD
RETURN
ENTRY TPHLDG(DF,DI,DA,DC,DD)
DF=DFST
DI=DMIN
DA=DMAX
DC=DFAC
DD=DVID
RETURN
END
C SUBROUTINE TPCSEC ====*====3====*====4====*====5====*====6====*====7
C
C EXTRA ENTRANCE FOR TPERSP
C WHEN CROS SECTION IS PLOT IN 2 DIMENSIONAL FIGURE
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPCSEC(XFU,XFG)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW,IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
DIMENSION XFU(3,2),XFG(3)
WRITE(6,*) ' ----- WELCOME TO TPERSP V3.0 2001/07/07 -----'
DO 10 I=1,3
IF(XFU(I,1).GE.XFU(I,2)) GO TO 999
IF(XFG(I).LE.0.0) GO TO 999
10 CONTINUE
IUOD=1
IHL=-2
IPR=0
IMOD=1
X0=XFU(1,1)
Y0=XFU(2,1)
Z0=XFU(3,1)
XN=XFU(1,2)
YN=XFU(2,2)
ZN=XFU(3,2)
XW=XFG(1)
YW=XFG(2)
ZW=XFG(3)
XFCTR=XW/(XN-X0)
YFCTR=YW/(YN-Y0)
ZFCTR=ZW/(ZN-Z0)
WRITE(6,600) X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW
WRITE(6,603) XFCTR,YFCTR,ZFCTR
600 FORMAT(' YOUR FIGURE IS FROM (',3E14.5,')'
& /' TO (',3E14.5,')'
& /' SIZE OF FIGURE (',3E14.5,')')
603 FORMAT(' SETTING FACTOR (',3E14.5,')')
RETURN
999 CONTINUE
WRITE(6,690)
690 FORMAT(' *** INVALID VALUE OF ARGUMENT IN TPCSEC ***')
RETURN
END
C SUBROUTINE TPLINE ====*====3====*====4====*====5====*====6====*====7
C
C THREE DIMENSIONAL LINE
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPLINE(X1,Y1,Z1,X2,Y2,Z2,NC)
IMPLICIT REAL*8(A-H,O-Z)
IF(NC.LT.1) RETURN
DX=(X2-X1)/FLOAT(NC)
DY=(Y2-Y1)/FLOAT(NC)
DZ=(Z2-Z1)/FLOAT(NC)
CALL TPPLOT(X1,Y1,Z1,3)
CALL TPPLOT(X1,Y1,Z1,2)
DO 1 I=1,NC
XX=X1+DX*FLOAT(I)
YY=Y1+DY*FLOAT(I)
ZZ=Z1+DZ*FLOAT(I)
CALL TPPLOT(XX,YY,ZZ,2)
1 CONTINUE
RETURN
END
C SUBROUTINE TPPOSI ====*====3====*====4====*====5====*====6====*====7
C
C SCREEN POINT XP, ZP CORRESPOND TO THREE DIMENSIONAL POINT
C XA,YA,ZA
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPPOSI(XA,YA,ZA,XP,ZP)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW,IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
IF(IMOD.EQ.1) GO TO 1
X=XA
Y=YA
Z=ZA
GO TO 4
1 X=(XA-X0)*XFCTR
Y=(YA-Y0)*YFCTR
Z=(ZA-Z0)*ZFCTR
4 CONTINUE
V=1.0/(X*SA+Y*CA+ET)
XP=(X*CA-Y*SA+SG)*V*EL-XP0
ZP=((Z+ZC)*V-ZM)*EL
RETURN
END
C SUBROUTINE TPXFXG ====*====3====*====4====*====5====*====6====*====7
C
C TRANSFORMATION FROM FUNCTION UNIT TO REAL UNIT
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPXFXG(XF,XG)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
REAL*8 XF(3),XG(3)
XG(1)=(XF(1)-XM(1,1))*XFCTR
XG(2)=(XF(2)-XM(2,1))*YFCTR
XG(3)=(XF(3)-XM(3,1))*ZFCTR
RETURN
END
C SUBROUTINE TPXGXF ====*====3====*====4====*====5====*====6====*====7
C
C TRANSFORMATION REAL UNIT TO FUNCTION UNIT
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPXGXF(XG,XF)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
REAL*8 XF(3),XG(3)
XF(1)=XG(1)/XFCTR+XM(1,1)
XF(2)=XG(2)/YFCTR+XM(2,1)
XF(3)=XG(3)/ZFCTR+XM(3,1)
RETURN
END
C SUBROUTINE TPPLOT ====*====3====*====4====*====5====*====6====*====7
C
C TRANSFORMATION FROM 3 DIMENSION TO CORDINATES ON THE SCREEN
C PEN UP AND DOWN IS CONTROLED BY TPHILN(ENTRY OF TPHILP)
C LINE FORM IS CONTROLED BY TPHILM(ENTRY OF TPHILD)
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPPLOT(XA,YA,ZA,JP)
C **** JP=3...START ****
C **** JP=0 OR 2...CONTINUED ****
IMPLICIT REAL*8(A-H,O-Z)
SAVE KPL,KPLD,IPC,XP,ZP
COMMON/TPERS1/X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW,IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
COMMON/TPERS8/X3,X4,Y3,Y4,Z3,Z4
& ,XPS1,ZPS1,XP1,ZP1
& ,X5,X6,Y5,Y6,Z5,Z6
& ,XPDS1,ZPDS1,XPDE1,ZPDE1
DIMENSION Q(3)
IF(IMOD.EQ.1) GO TO 1
X=XA
Y=YA
Z=ZA
XX=(X/XFCTR)+X0
YY=(Y/YFCTR)+Y0
ZZ=(Z/ZFCTR)+Z0
GO TO 4
1 XX=XA
YY=YA
ZZ=ZA
X=(XX-X0)*XFCTR
Y=(YY-Y0)*YFCTR
Z=(ZZ-Z0)*ZFCTR
4 CONTINUE
IM=IMOD
IMOD=-1
IP=IPC
XPS1=XP
ZPS1=ZP
KPLS=KPL
KPLDS=KPLD
V=1.0/(X*SA+Y*CA+ET)
XP=(X*CA-Y*SA+SG)*V*EL-XP0
ZP=((Z+ZC)*V-ZM)*EL
XQ=X
YQ=Y
ZQ=Z
Q(1)=XQ
Q(2)=YQ
Q(3)=ZQ
IF(JP.NE.3) GO TO 2
C **** JP=3 ****
CALL TPHILN(Q,KPL)
CALL TPHILM(Q,KPLD,IS1,IS2,IS3)
IF(KPL.EQ.-1) IPC=3
IF(KPL.EQ.1) IPC=2
IF(KPLD.EQ.-1) CALL SETDAS(IS1,IS2,IS3)
IF(KPLD.EQ.1) CALL SETDAS(0,0,0)
IF(IPR.EQ.2) WRITE(6,600) XX,YY,ZZ,X,Y,Z,XP,ZP,JP,IPC,KPL,KPLD
600 FORMAT(' TRACE IN TPPLOT',3E11.3,5F10.4,4I4)
CALL PLOTPS(XP,ZP,2)
GO TO 800
C **** JP=0 OR 2 ****
2 CONTINUE
CALL TPHILN(Q,KPL)
CALL TPHILM(Q,KPLD,IS1,IS2,IS3)
IF(KPL.EQ.KPLS.AND.KPLD.EQ.KPLDS) GO TO 3
X3=X4
X4=X
Y3=Y4
Y4=Y
Z3=Z4
Z4=Z
XP1=XP
ZP1=ZP
CALL TPVANS(IP,KPLS,KPL,KPLDS,KPLD)
IPC=IP
GO TO 800
C ****
3 CONTINUE
IPV=IP
IF(IPR.EQ.2) WRITE(6,600) XX,YY,ZZ,X,Y,Z,XP,ZP,JP,IP,KPL,KPLD
CALL PLOTPS(XP,ZP,IPV-1)
800 CONTINUE
X4=X
Y4=Y
Z4=Z
IMOD=IM
RETURN
END
C SUBROUTINE TPVANS ====*====3====*====4====*====5====*====6====*====7
C
C FIND THE VANISHING POINT
C PEN UP AND DOWN IS CONTROLED BY TPHILN
C LINE FORM IS CONTROLED BY TPHILM
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPVANS(KP,KPLS,KPLE,KPLDS,KPLDE)
C ** TO FIND VANISH POINT **
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/W(9),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
COMMON/TPERS8/X3,X4,Y3,Y4,Z3,Z4,
& XPS1,ZPS1,XP1,ZP1
& ,X5,X6,Y5,Y6,Z5,Z6
& ,XPDS1,ZPDS1,XPDE1,ZPDE1
DIMENSION Q(3)
DATA PX,PZ,PX1,PZ1/0.05,0.05,0.055,0.055/
IF(KPLE.EQ.KPLS) THEN
X5=X3
Y5=Y3
Z5=Z3
X6=X4
Y6=Y4
Z6=Z4
XPDS1=XPS1
ZPDS1=ZPS1
XPDE1=XP1
ZPDE1=ZP1
CALL TPVAND(KP,KPLS,KPLDE)
RETURN
END IF
DX=-0.5*(X4-X3)
DY=-0.5*(Y4-Y3)
DZ=-0.5*(Z4-Z3)
XP=XP1
ZP=ZP1
KPL=KPLE
KPL1=KPLE
XQ=X4
YQ=Y4
ZQ=Z4
DO 130 L=1,8
XQ=XQ+DX
YQ=YQ+DY
ZQ=ZQ+DZ
V=1.0/(XQ*SA+YQ*CA+ET)
XPS=XP
ZPS=ZP
XP=(XQ*CA-YQ*SA+SG)*V*EL-XP0
ZP=((ZQ+ZC)*V-ZM)*EL
KPLS=KPL
Q(1)=XQ
Q(2)=YQ
Q(3)=ZQ
CALL TPHILN(Q,KPL)
IF(ABS(XP-XPS).GT.PX) GO TO 131
IF(ABS(ZP-ZPS).GT.PZ) GO TO 131
GO TO 132
131 CONTINUE
DX=0.5*DX
DY=0.5*DY
DZ=0.5*DZ
IF(KPL.EQ.KPLS) GO TO 130
DX=-DX
DY=-DY
DZ=-DZ
130 CONTINUE
132 CONTINUE
IF(ABS(XP-XP1).GT.PX1)
& GO TO 133
IF(ABS(ZP-ZP1).GT.PZ1)
& GO TO 133
XP=XP1
ZP=ZP1
XQ=X4
YQ=Y4
ZQ=Z4
GO TO 134
133 CONTINUE
IF(ABS(XP-XPS1).GT.PX1)
& GO TO 134
IF(ABS(ZP-ZPS1).GT.PZ1)
& GO TO 134
XP=XPS1
ZP=ZPS1
XQ=X3
YQ=Y3
ZQ=Z3
134 CONTINUE
Q(1)=XQ
Q(2)=YQ
Q(3)=ZQ
CALL TPHILM(Q,KPLD,IS1,IS2,IS3)
IF(KPLD.NE.KPLDS) THEN
X5=X3
Y5=Y3
Z5=Z3
X6=XQ
Y6=YQ
Z6=ZQ
XPDS1=XPS1
ZPDS1=ZPS1
XPDE1=XP
ZPDE1=ZP
CALL TPVAND(KP,KPLS,KPLD)
ELSE
KPV=KP
IF(IPR.EQ.2) WRITE(6,600) XQ,YQ,ZQ,XP,ZP,KP,KPL,KPLD
600 FORMAT(' TRACE IN TPVANS',33X,5F10.4,4X,3I4)
CALL PLOTPS(XP,ZP,KPV-1)
END IF
KP=5-KP
IF(KPLD.NE.KPLDE) THEN
X5=XQ
Y5=YQ
Z5=ZQ
X6=X4
Y6=Y4
Z6=Z4
XPDS1=XP
ZPDS1=ZP
XPDE1=XP1
ZPDE1=ZP1
CALL TPVAND(KP,KPLE,KPLDE)
ELSE
KPL=KPL1
XP=XP1
ZP=ZP1
KPV=KP
IF(IPR.EQ.2) WRITE(6,600) X4,Y4,Z4,XP,ZP,KP,KPL,KPLDE
CALL PLOTPS(XP,ZP,KPV-1)
END IF
C ****
RETURN
END
C SUBROUTINE TPVAND ====*====3====*====4====*====5====*====6====*====7
C
C FIND THE CHANGE POINT OF LINE FORM
C LINE FORM IS CONTROLED BY TPHILM
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPVAND(KP,KPLH,KPLDE)
C ** TO FIND CHANGE POINT **
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/W(9),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
COMMON/TPERS8/X3,X4,Y3,Y4,Z3,Z4,
& XPS1,ZPS1,XP1,ZP1
& ,X5,X6,Y5,Y6,Z5,Z6
& ,XPDS1,ZPDS1,XPDE1,ZPDE1
DIMENSION Q(3)
DATA PX,PZ,PX1,PZ1/0.05,0.05,0.055,0.055/
DX=-0.5*(X6-X5)
DY=-0.5*(Y6-Y5)
DZ=-0.5*(Z6-Z5)
XP=XPDE1
ZP=ZPDE1
KPL1=KPLDE
KPL=KPLDE
XQ=X6
YQ=Y6
ZQ=Z6
DO 130 L=1,8
XQ=XQ+DX
YQ=YQ+DY
ZQ=ZQ+DZ
V=1.0/(XQ*SA+YQ*CA+ET)
XPS=XP
ZPS=ZP
XP=(XQ*CA-YQ*SA+SG)*V*EL-XP0
ZP=((ZQ+ZC)*V-ZM)*EL
KPLS=KPL
Q(1)=XQ
Q(2)=YQ
Q(3)=ZQ
CALL TPHILM(Q,KPL,IS1,IS2,IS3)
IF(ABS(XP-XPS).GT.PX) GO TO 131
IF(ABS(ZP-ZPS).GT.PZ) GO TO 131
GO TO 132
131 CONTINUE
DX=0.5*DX
DY=0.5*DY
DZ=0.5*DZ
IF(KPL.EQ.KPLS) GO TO 130
DX=-DX
DY=-DY
DZ=-DZ
130 CONTINUE
132 CONTINUE
IF(ABS(XP-XPDE1).GT.PX1)
& GO TO 133
IF(ABS(ZP-ZPDE1).GT.PZ1)
& GO TO 133
XP=XPDE1
ZP=ZPDE1
XQ=X6
YQ=Y6
ZQ=Z6
GO TO 134
133 CONTINUE
IF(ABS(XP-XPDS1).GT.PX1)
& GO TO 134
IF(ABS(ZP-ZPDS1).GT.PZ1)
& GO TO 134
XP=XPDS1
ZP=ZPDS1
XQ=X5
YQ=Y5
ZQ=Z5
134 CONTINUE
KPV=KP
IF(IPR.EQ.2) WRITE(6,600) XQ,YQ,ZQ,XP,ZP,KP,KPLH,KPL
600 FORMAT(' TRACE IN TPVAND',33X,5F10.4,4X,3I4)
CALL PLOTPS(XP,ZP,KPV-1)
IF(KPLDE.EQ.-1) CALL SETDAS(IS1,IS2,IS3)
C WRITE(6,*) ' SS1 SS2 NLFF ',SS1,SS2,NLFF
IF(KPLDE.EQ.1) CALL SETDAS(0,0,0)
XP=XPDE1
ZP=ZPDE1
KPLD=KPLDE
KPV=KP
IF(IPR.EQ.2) WRITE(6,600) X4,Y4,Z4,XP,ZP,KP,KPLH,KPLD
CALL PLOTPS(XP,ZP,2)
CALL PLOTPS(XP,ZP,KPV-1)
C ****
RETURN
END
C SUBROUTINE TPHILP ====*====3====*====4====*====5====*====6====*====7
C
C FIRST HIDDEN LINE ROUTINE
C IHL=0 -------- NORMAL CASE
C IF(THERE IS ANY POINT(XA) ON THE WAY FROM Q TO EYE POINT
C FUNC(XA)-HHV)*IUOD.LT.0.0)
C OR POINT Q IS OUTSIDE
C OR POINT Q IS CUT OUT BY TPHILL THEN KPL=-1
C ELSE KPL=1
C IHL=2
C IF POINT Q IS OUTSIDE
C OR POINT Q IS CUT OUT BY TPHILL THEN KPL=-1
C ELSE KPL=1
C IHL=3
C ALMOST THE SAME AS IHL=0
C BUT EVEN IF(POINT Q IS CUT OUT BY TPHILL)
C KPL IS NOT AUTOMATICALLY -1
C IHL=4, KPL IS ALWAYS 1,UNLESS Q IS OUTtSIDE
C IHL=5 KPL IS ALWAYS 1
C INITIAL SET IN TPERSP
C IIOD=1
C IHL =0
C CALL TPHILP(-1.0)
C 1989/12/09 BY A.YANASE
C 2001/03/09/ BY A.YANASE
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPHILP(HHV)
IMPLICIT REAL*8(A-H,O-Z)
EXTERNAL FUNC,LFUNC
COMMON/TPERS1/X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW,IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,W(8),XEYE,YEYE,ZEYE
COMMON/TPERS6/DFST,DMIN,DMAX,DFAC,DVID
DIMENSION Q(3),XA(3),P(3),QQ(3)
HV=HHV
RETURN
C
ENTRY TPHILN(Q,KPL)
KPL=1
IF(IHL.EQ.5) GO TO 2
IF(IMOD.EQ.-1) GO TO 11
XX=Q(1)
YY=Q(2)
ZZ=Q(3)
XQ=(XX-X0)*XFCTR
YQ=(YY-Y0)*YFCTR
ZQ=(ZZ-Z0)*ZFCTR
GO TO 12
11 XQ=Q(1)
YQ=Q(2)
ZQ=Q(3)
XX=(XQ/XFCTR)+X0
YY=(YQ/YFCTR)+Y0
ZZ=(ZQ/ZFCTR)+Z0
12 CONTINUE
IF(XQ.LT.-1.0D-3) GO TO 4
IF(XQ.GT.XW+1.0D-3) GO TO 4
IF(YQ.LT.-1.0D-3) GO TO 4
IF(YQ.GT.YW+1.0D-3) GO TO 4
IF(ZQ.LT.-1.0D-3) GO TO 4
IF(ZQ.GT.ZW+1.0D-3) GO TO 4
IF(IHL.EQ.4) GO TO 2
XA(1)=XX
XA(2)=YY
XA(3)=ZZ
LL=LFUNC(XA)
C WRITE(6,600) LL,XX,YY,ZZ,XQ,YQ,ZQ
C 600 FORMAT(I4,6E13.5)
IF(LL.EQ.1.AND.IHL.NE.3) GO TO 4
IF(IHL.EQ.2) GO TO 2
7 CONTINUE
PX=XEYE-XQ
PY=YEYE-YQ
PZ=ZEYE-ZQ
WA=SQRT(PX*PX+PY*PY+PZ*PZ)
PX=PX/WA
PY=PY/WA
PZ=PZ/WA
IF(IHL.EQ.1) GO TO 123
DD=DFST
DST=DMIN
DW=1.0D30
DO 1 I=1,20000
XD=XQ+PX*DD
IF(XD.LT.0.0) GO TO 2
IF(XD.GT.XW) GO TO 2
YD=YQ+PY*DD
IF(YD.LT.0.0) GO TO 2
IF(YD.GT.YW) GO TO 2
ZD=ZQ+PZ*DD
IF(ZD.LT.0.0) GO TO 2
IF(ZD.GT.ZW) GO TO 2
XA(1)=(XD/XFCTR)+X0
XA(2)=(YD/YFCTR)+Y0
XA(3)=(ZD/ZFCTR)+Z0
IF(LFUNC(XA).EQ.1) GO TO 6
EE=FUNC(XA)
C *******************************
C WRITE(6,601) EE,HV,XX,YY,ZZ,XD,YD,ZD,DD
C 601 FORMAT(9E12.4)
C *****************************
IF(EE.GT.HV) GO TO 3
IF(IUOD.EQ.-1) GO TO 5
GO TO 4
3 IF(IUOD.EQ.1) GO TO 5
GO TO 4
5 DV=ABS(EE-HV)
C *******************************************************
IF(DV.LT.DW.AND.DST*DVID.GE.DMIN) DST=DST*DVID
C IF(DV.GT.DW.AND.DST*DFAC.LE.DMAX) DST=DST*DFAC
IF(DV.GE.DW.AND.DST*DFAC.LE.DMAX) DST=DST*DFAC
C **********************************************************
DW=DV
6 DD=DD+DST
1 CONTINUE
4 KPL=-1
2 CONTINUE
RETURN
123 CONTINUE
P(1)=PX
P(2)=PY
P(3)=PZ
QQ(1)=XX
QQ(2)=YY
QQ(3)=ZZ
EE=TPMASB(QQ,P)
IF(EE.LE.0.0) KPL=-1
RETURN
END
C SUBROUTINE TPHILD ====*====3====*====4====*====5====*====6====*====7
C
C SECOND HIDDEN LINE ROUTINE
C IF(THERE IS ANY POINT(XA) ON THE WAY FROM Q TO EYE POINT
C FUNC(XA)-HHV)*IUOD.LT.0.0)
C OR POINT Q IS OUTSIDE
C OR POINT Q IS CUT OUT BY TPHILL THEN KPL=-1
C ELSE KPL=1
C HOWEVER IF(IFLAG.EQ.0) KPL IS ALWAYS -1
C INITIAL SET IN TPERSP
C IUOD=1
C CALL TPHILD(0.0,0,0,0,0)
C PART OF KPL=-1 IS PLOTTED BY THE LINE GIVEN BY
C IS1,IS2,IS3
C IS1=0 FULL LINE
C IS3=0 DASHED LINE
C IS3.GT.0 DASH AND DOT LINE
C ---------- ---- -----------
C IS1 IS2 IS3 IS2 IS1
C 1989/12/09 BY A.YANASE
C 2001/03/09 BY A.YANASE
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPHILD(HHV,IFLAG,IS1,IS2,IS3)
IMPLICIT REAL*8(A-H,O-Z)
SAVE ISS1,ISS2,ISS3,HV,IIFL
COMMON/TPERS1/X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW,IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,W(8),XEYE,YEYE,ZEYE
COMMON/TPERS6D/DFST,DMIN,DMAX,DFAC,DVID
DIMENSION Q(3),XA(3)
HV=HHV
IIFL=IFLAG
ISS1=IS1
ISS2=IS2
ISS3=IS3
RETURN
C
C
ENTRY TPHILM(Q,KPL,IS1,IS2,IS3)
IS1=ISS1
IS2=ISS2
IS3=ISS3
IF(IIFL.EQ.0) THEN
KPL=-1
RETURN
END IF
KPL=1
IF(IMOD.EQ.-1) GO TO 11
XX=Q(1)
YY=Q(2)
ZZ=Q(3)
XQ=(XX-X0)*XFCTR
YQ=(YY-Y0)*YFCTR
ZQ=(ZZ-Z0)*ZFCTR
GO TO 12
11 XQ=Q(1)
YQ=Q(2)
ZQ=Q(3)
XX=(XQ/XFCTR)+X0
YY=(YQ/YFCTR)+Y0
ZZ=(ZQ/ZFCTR)+Z0
12 CONTINUE
XA(1)=XX
XA(2)=YY
XA(3)=ZZ
PX=XEYE-XQ
PY=YEYE-YQ
PZ=ZEYE-ZQ
WA=SQRT(PX*PX+PY*PY+PZ*PZ)
PX=PX/WA
PY=PY/WA
PZ=PZ/WA
DD=DFST
DST=DMIN
DW=1.0E30
DO 1 I=1,20000
XD=XQ+PX*DD
IF(XD.LT.0.0) GO TO 2
IF(XD.GT.XW) GO TO 2
YD=YQ+PY*DD
IF(YD.LT.0.0) GO TO 2
IF(YD.GT.YW) GO TO 2
ZD=ZQ+PZ*DD
IF(ZD.LT.0.0) GO TO 2
IF(ZD.GT.ZW) GO TO 2
XA(1)=(XD/XFCTR)+X0
XA(2)=(YD/YFCTR)+Y0
XA(3)=(ZD/ZFCTR)+Z0
IF(LFUNC(XA).EQ.1) GO TO 6
EE=FUNCD(XA)
C *******************************
C WRITE(6,601) EE,HV,XX,YY,ZZ,XD,YD,ZD,DD
C 601 FORMAT(9E12.4)
C *****************************
IF(EE.GT.HV) GO TO 3
IF(IUOD.EQ.-1) GO TO 5
GO TO 4
3 IF(IUOD.EQ.1) GO TO 5
GO TO 4
5 DV=ABS(EE-HV)
C *******************************************************
IF(DV.LT.DW.AND.DST*DVID.GE.DMIN) DST=DST*DVID
C IF(DV.GT.DW.AND.DST*DFAC.LE.DMAX) DST=DST*DFAC
IF(DV.GE.DW.AND.DST*DFAC.LE.DMAX) DST=DST*DFAC
C **********************************************************
DW=DV
6 DD=DD+DST
1 CONTINUE
4 KPL=-1
2 CONTINUE
RETURN
END
C SUBROUTINE TPSECT ====*====3====*====4====*====5====*====6====*====7
C
C THREE DIMENSIONAL(IHL.GE.0) OR
C TWO DIMENSIONAL(IHL.LT.0) PLOT OF
C CROSS SECTION OF AA=FUNC(X,Y,Z) BY THE PLAN
C THROUGH XCC WITH DIRCTION COSINE GIVEN BY DRR
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPSECT(FUNC,AA,DRR,XCC,IAL,ICR,NNX,NNY,KKX,KKY)
IMPLICIT REAL*8(A-H,O-Z)
EXTERNAL FUNC
REAL*8 KB,KBM
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W2(11)
COMMON/TPERS3/ALT(2402,2402),HX,HY,NX,NY,HV
COMMON/TPERS4/XOR(3),WW(3,3),XFC,YFC
DIMENSION DR(3),XCD(3),XC(3),DRR(3),XCC(3)
DIMENSION XYZJ(3),XYZI(3)
DIMENSION DW(3),XMAX(2),XMIN(2),XCR(3,12)
DIMENSION KB(3),KBM(3,4),XD(3),WX(3)
IF(NNX.GT.2402) GO TO 999
IF(NNX.LE.0) GO TO 999
IF(NNY.GT.2402) GO TO 999
IF(NNY.LE.0) GO TO 999
CALL PTIME(TB)
NX=NNX
NY=NNY
HV=AA
IM=IMOD
IMOD=-1
DO 11 I=1,3
DR(I)=DRR(I)
11 XC(I)=XCC(I)
IF(IM.EQ.1) THEN
DO 12 I=1,3
DR(I)=DR(I)*FCTR(I)
12 XC(I)=(XC(I)-XM(I,1))*FCTR(I)
END IF
IF(IAL.EQ.1) GO TO 810
WA=DR(1)*DR(1)+DR(2)*DR(2)+DR(3)*DR(3)
WA=SQRT(WA)
DL=DR(1)/WA
DM=DR(2)/WA
DN=DR(3)/WA
DW(1)=DL
DW(2)=DM
DW(3)=DN
WA=SQRT(DL*DL+DM*DM)
IF(WA.LT.1.0D-10) THEN
DO 18 I=1,3
DO 19 J=1,3
19 WW(I,J)=0.0
18 WW(I,I)=1.0
ELSE
WW(1,1)=DM/WA
WW(1,2)=(DN*DL)/WA
WW(1,3)=DL
WW(2,1)=-DL/WA
WW(2,2)=(DN*DM)/WA
WW(2,3)=DM
WW(3,1)=0.0
WW(3,2)=-WA
WW(3,3)=DN
END IF
C WRITE(6,602) XM,WW,XC,DW
C 602 FORMAT(9E14.6)
WA=DW(1)*XC(1)+DW(2)*XC(2)+DW(3)*XC(3)
C EQUATION OF PLANE IS GIVEN BY
C WA=DW(1)*X+DW(2)*Y+DW(3)*Z
NCR=0
XMAX(1)=-1.0E20
XMIN(1)=1.0E20
XMAX(2)=-1.0E20
XMIN(2)=1.0E20
DO 20 I=1,2
DO 21 J=I+1,3
K=J-I
IF(J.EQ.2) K=3
C I=1, J=2, K=3
C I=1, J=3, K=2
C I=2, J=3, K=1
C IF(IPR.GE.1) WRITE(6,600) I,J,K
IF(ABS(DW(K)).LT.1.0D-10) GO TO 21
DO 22 I1=1,2
IF(I1.EQ.1) THEN
XYZI(I)=0.0
XYZI(J)=0.0
XYZI(K)=0.0
ELSE
XYZI(I)=XM(I,3)
XYZI(J)=XM(J,3)
XYZI(K)=XM(K,3)
END IF
DO 23 J1=1,2
IF(J1.EQ.1) THEN
XYZJ(I)=0.0
XYZJ(J)=0.0
XYZJ(K)=0.0
ELSE
XYZJ(I)=XM(I,3)
XYZJ(J)=XM(J,3)
XYZJ(K)=XM(K,3)
END IF
XH=(WA-DW(I)*XYZI(I)-DW(J)*XYZJ(J))/DW(K)
IF(XH.LT.-1.0D-4) GO TO 23
IF(XH.GT.XM(K,3)+1.0D-4) GO TO 23
WX(I)=XYZI(I)
WX(J)=XYZJ(J)
WX(K)=XH
NCR=NCR+1
C WRITE(6,600) I,J,K,WX
C 600 FORMAT(3I10,3E20.7)
DO 24 II=1,3
WB=0.0
DO 25 JJ=1,3
25 WB=WB+WX(JJ)*WW(JJ,II)
XCR(II,NCR)=WB
IF(WB.GT.XMAX(II).AND.II.NE.3) XMAX(II)=WB
IF(WB.LT.XMIN(II).AND.II.NE.3) XMIN(II)=WB
24 CONTINUE
23 CONTINUE
22 CONTINUE
21 CONTINUE
20 CONTINUE
IF(XMAX(1).LE.XMIN(1)) GO TO 999
IF(XMAX(2).LE.XMIN(2)) GO TO 999
C WRITE(6,600) NCR
C WRITE(6,602) XCR
C WRITE(6,601) XMAX,XMIN
C 601 FORMAT(4E20.7)
XCD(1)=(XMAX(1)+XMIN(1))*0.5
XCD(2)=(XMAX(2)+XMIN(2))*0.5
XCD(3)=XCR(3,1)
C WRITE(6,602) XCD
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6001) TE,XC,DW
6001 FORMAT(' TIME FOR SETTING TPSECT',F7.4
& /10X,' XC',3F9.5,' DR',3F9.5)
TB=TA
HX=(XMAX(1)-XMIN(1))/FLOAT(NX-1)
HY=(XMAX(2)-XMIN(2))/FLOAT(NY-1)
NXH=(NX-1)/2+1
NYH=(NY-1)/2+1
HMAX=-1.0E20
HMIN=1.0E20
XD(3)=XCD(3)
DO 13 I=1,NX
XD(1)=XCD(1)+FLOAT(I-NXH)*HX
DO 14 J=1,NY
XD(2)=XCD(2)+FLOAT(J-NYH)*HY
DO 15 I1=1,3
WA=0.0
DO 16 I2=1,3
16 WA=WA+WW(I1,I2)*XD(I2)
KB(I1)=WA/FCTR(I1)+XM(I1,1)
15 CONTINUE
WB=FUNC(KB)
ALT(I,J)=WB
IF(WB.GT.HMAX) HMAX=WB
IF(WB.LT.HMIN) HMIN=WB
IF(J.EQ.1) THEN
IF(I.EQ.1) THEN
DO 33 I1=1,3
XOR(I1)=(KB(I1)-XM(I1,1))*FCTR(I1)
33 KBM(I1,1)=KB(I1)
ELSE IF(I.EQ.NX) THEN
DO 32 I1=1,3
32 KBM(I1,2)=KB(I1)
END IF
ELSE IF(J.EQ.NY) THEN
IF(I.EQ.1) THEN
DO 31 I1=1,3
31 KBM(I1,3)=KB(I1)
ELSE IF(I.EQ.NX) THEN
DO 30 I1=1,3
30 KBM(I1,4)=KB(I1)
END IF
END IF
14 CONTINUE
13 CONTINUE
IF(IPR.GE.1) THEN
WRITE(6,6020)
6020 FORMAT(' CORNER COORDINATES IN FUNCTION UNIT OF SETTING PLANE')
DO 29 J1=1,4
WRITE(6,6021) (KBM(I1,J1),I1=1,3)
6021 FORMAT(1H ,3F10.5)
29 CONTINUE
END IF
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6002) TE,HMAX,HV,HMIN
6002 FORMAT(' TIME FOR VALUES TPSECT ',F7.4
& ,3E11.4)
TB=TA
IF(HMAX.LT.HV) GO TO 820
IF(HMIN.GT.HV) GO TO 820
IF(IHL.GE.0.AND.ICR.NE.1) GO TO 810
WB=0.0
WA=0.0
DO 34 I1=1,3
WC=(KBM(I1,1)-KBM(I1,2))**2
WA=WA+WC
WB=WB+WC*(FCTR(I1)**2)
34 CONTINUE
XFC=SQRT(WB/WA)
WWX=SQRT(WB)
WB=0.0
WA=0.0
DO 35 I1=1,3
WC=(KBM(I1,1)-KBM(I1,3))**2
WA=WA+WC
WB=WB+WC*(FCTR(I1)**2)
35 CONTINUE
YFC=SQRT(WB/WA)
WWY=SQRT(WB)
IF(IHL.GE.0) GO TO 810
CALL PLOTPS(0.0D0,0.0D0,2)
CALL PLOTPS(WWX,0.0D0,1)
CALL PLOTPS(WWX,WWY,1)
CALL PLOTPS(0.0D0,WWY,1)
CALL PLOTPS(0.0D0,0.0D0,1)
810 CONTINUE
CALL TPCONT
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6003) TE
6003 FORMAT(' CPTIME FOR TPCONT IN TPSECT',F7.4)
TB=TA
IF(ICR.NE.1) GO TO 820
NCX=KKX
NCY=KKY
CALL TPCROS(NCX,NCY)
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6004) TE
6004 FORMAT(' CPTIME FOR TPCROS IN TPSECT',F7.4)
TB=TA
820 IMOD=IM
RETURN
999 CONTINUE
WRITE(6,690)
690 FORMAT(' *** INVALID VALUE OF ARGUMENT IN TPSECT ***')
IMOD=IM
RETURN
END
C SUBROUTINE TPFXYZ ====*====3====*====4====*====5====*====6====*====7
C
C THREE DIMENSIONAL PLOT OF AA=FUNC(X,Y,Z)
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPFXYZ(FUNC,AA,NN,NK,ICRO)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 KB,KBA
EXTERNAL FUNC
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W2(11)
COMMON/TPERS3/ALT(2402,2402),HX,HY,NX,NY,EF
COMMON/TPERS4/XOR(3),WW(9),XFC,YFC
DIMENSION KB(3),KBA(3),NN(3),NK(3)
DIMENSION W(9,3)
DATA W/1.0,0.0,0.0, 0.0,1.0,0.0, 0.0,0.0,1.0,
& 0.0,1.0,0.0, 0.0,0.0,1.0, 1.0,0.0,0.0,
& 0.0,0.0,1.0, 1.0,0.0,0.0, 0.0,1.0,0.0/
DO 51 I=1,3
IF(NN(I).GT.2402) GO TO 999
IF(NN(I).LE.0) GO TO 999
IF(NK(I).LE.0) GO TO 999
51 CONTINUE
IM=IMOD
IMOD=1
CALL PTIME(TB)
EF=AA
ICROS=ICRO
DO 14 IXX=1,3
IYY=MOD(IXX,3)+1
IZZ=MOD(IXX+1,3)+1
HX=(XM(IXX,2)-XM(IXX,1))/FLOAT(NN(IXX)-1)
HY=(XM(IYY,2)-XM(IYY,1))/FLOAT(NN(IYY)-1)
HZ=(XM(IZZ,2)-XM(IZZ,1))/FLOAT(NN(IZZ)-1)
DO 63 I=1,9
WW(I)=W(I,IXX)
63 CONTINUE
XOR(IXX)=XM(IXX,1)
XOR(IYY)=XM(IYY,1)
NX=NN(IXX)
NY=NN(IYY)
NZ=NN(IZZ)
NKZ=NK(IZZ)
NCX=NK(IXX)
NCY=NK(IYY)
DO 11 KZ=1,NZ,NKZ
KB(IZZ)=XM(IZZ,1)+HZ*FLOAT(KZ-1)
XOR(IZZ)=KB(IZZ)
EMAX=-1.0E38
EMIN=1.0E38
DO 21 KY=1,NY
KB(IYY)=XM(IYY,1)+HY*FLOAT(KY-1)
DO 22 KX=1,NX
KB(IXX)=XM(IXX,1)+HX*FLOAT(KX-1)
DO 12 I=1,3
12 KBA(I)=KB(I)
WA=FUNC(KBA)
ALT(KX,KY)=WA
IF(EMAX.LT.WA) EMAX=WA
IF(EMIN.GT.WA) EMIN=WA
22 CONTINUE
21 CONTINUE
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6003) TE,IXX,IYY,IZZ,KB(IZZ)
& ,EMAX,EF,EMIN
6003 FORMAT(' TIME FOR VALUES TPFXYZ',F7.4,3I2,4E11.4)
TB=TA
IF(EMAX.LT.EF) GO TO 11
IF(EMIN.GT.EF) GO TO 11
CALL TPCONT
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6001) TE
6001 FORMAT(' TIME FOR TPCONT TPFXYZ',F7.3)
TB=TA
IF(ICROS.NE.1) GO TO 11
IF(KZ.NE.1.AND.KZ.NE.NZ) GO TO 11
CALL TPCROS(NCX,NCY)
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6002) TE
6002 FORMAT(' TIME FOR TPCROS TPFXYZ',F7.4)
TB=TA
11 CONTINUE
14 CONTINUE
IMOD=IM
RETURN
999 CONTINUE
WRITE(6,690)
690 FORMAT(' *** INVALID VALUE OF ARGUMENT IN TPFXYZ ***')
RETURN
END
C SUBROUTINE TPCONT ====*====3====*====4====*====5====*====6====*====7
C
C TPCONT
C
C 1989/12/14 A. YANASE
C ORIGINALY BY S.MIYAHARA FOR AREA CALCULATION
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPCONT
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS3/ALT(2402,2402),HX,HY,NX,NY,HV
COMMON/TPERS5/KSW(2402,2402),AX,AY,BX,BY,
& UA1,UA2,KV,IREP,IA,JA,IT,JT,KS
C ****
IG=1
DO 50 I=1,NX
DO 50 J=1,NY
50 KSW(I,J)=0
C ****
KERR=0
KV=1
C ****
DO 20 KX=1,NX
V=ALT(KX,1)-HV
IF(V) 2,3,1
1 KG=1
GO TO 4
2 KG=2
GO TO 4
3 KG=3
4 CONTINUE
IF(KSW(KX,1).NE.KV) GO TO 5
IF(ALT(KX,2)-HV.LE.0.0) KG=2
5 CONTINUE
C ****
DO 30 KY=2,NY
V=ALT(KX,KY)-HV
IF(KSW(KX,KY).EQ.KV) GO TO 15
C ****
GO TO (6,8,21),KG
C ** ALT(KX,1),GT.HV **
6 IF(V) 7,7,30
7 CONTINUE
IREP=0
DO 60 M=1,2
IA=KX
JA=KY-1
GO TO (31,32),M
31 IF(KX.EQ.NX) GO TO 60
IT=1
KS=-1
GO TO 33
32 IF(KX.EQ.1) GO TO 60
IF(IREP.GE.2) GO TO 60
IT=-1
KS=1
KSW(IA,JA)=0
33 JT=0
UA1=ALT(IA,JA)
UA2=ALT(IA,JA+1)
ALX=HX*FLOAT(IA-1)
RT=(UA2-HV)/(UA1-UA2)
ALY=HY*(FLOAT(JA)+RT)
KERR=1
CALL TP2TO3(ALX,ALY,3)
CALL TP2TO3(ALX,ALY,2)
AX=ALX
AY=ALY
CALL TPSEAR
IF(KX.EQ.1.OR.KX.EQ.NX) GO TO 60
IF(ABS(AX-BX).GT.HX*2.2) GO TO 60
IF(ABS(AY-BY).GT.HY*2.2) GO TO 60
CALL TP2TO3(AX,AY,2)
IG=IG+1
GO TO 61
60 CONTINUE
61 CONTINUE
KG=2
GO TO 30
C ** ALT(KX,1).LT.HV **
8 IF(V) 30,30,9
9 CONTINUE
IREP=0
DO 70 M=1,2
IA=KX
JA=KY
GO TO (34,35),M
34 IF(KX.EQ.NX) GO TO 70
IT=1
KS=1
GO TO 36
35 IF(KX.EQ.1) GO TO 70
IF(IREP.GE.2) GO TO 70
IT=-1
KS=-1
KSW(IA,JA)=0
36 JT=0
UA1=ALT(IA,JA)
UA2=ALT(IA,JA-1)
ALX=HX*FLOAT(IA-1)
RT=(UA1-HV)/(UA1-UA2)
ALY=HY*(FLOAT(JA-1)-RT)
KERR=1
CALL TP2TO3(ALX,ALY,3)
CALL TP2TO3(ALX,ALY,2)
AX=ALX
AY=ALY
CALL TPSEAR
IF(KX.EQ.1.OR.KX.EQ.NX) GO TO 70
IF(ABS(AX-BX).GT.HX*2.2) GO TO 70
IF(ABS(AY-BY).GT.HY*2.2) GO TO 70
CALL TP2TO3(AX,AY,2)
GO TO 71
70 CONTINUE
71 CONTINUE
KG=1
GO TO 30
C ****
15 GO TO (16,17,30),KG
16 KG=2
IF(KY.EQ.NY) GO TO 30
IF(ALT(KX,KY+1)-HV.GT.0.0) KG=1
GO TO 30
17 KG=1
IF(V.LT.0.0) KG=2
GO TO 30
C ****
21 IF(V) 23,24,22
22 KG=1
GO TO 25
23 KG=2
GO TO 25
24 KG=3
25 CONTINUE
C ****
30 CONTINUE
C ****
20 CONTINUE
C ****
IF(KERR.EQ.0) WRITE(6,2001) HV
2001 FORMAT(' (TPAREB) HEIT (',E9.2,
& ') IS NOT FOUND.')
C ****
RETURN
END
C SUBROUTINE TPSEAR ====*====3====*====4====*====5====*====6====*====7
C
C TPSEAR
C
C 1989/12/14 A. YANASE
C ORIGINALY BY S.MIYAHARA FOR AREA CALCULATION
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPSEAR
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS3/ALT(2402,2402),HX,HY,NX,NY,HV
COMMON/TPERS5/KSW(2402,2402),AX,AY,BX,BY,
& UA1,UA2,KV,IREP,IA,JA,IT,JT,KS
C ****
ITV=0
JTV=0
LN=800
IREP=0
C ****
DO 10 L=2,LN
IF(KSW(IA,JA).EQ.KV) GO TO 806
GO TO 819
806 IF(IA+ITV.EQ.0.OR.JA+JTV.EQ.0) GO TO 99
IF(IA-ITV.EQ.0.OR.JA-JTV.EQ.0) GO TO 99
ALT1=ALT(IA-ITV,JA-JTV)
ALT2=ALT(IA,JA)
ALT3=ALT(IA+ITV,JA+JTV)
ALT4=ALT(IA-ITV,JA+JTV)
ALT5=ALT(IA+ITV,JA-JTV)
IF((ALT1.LT.ALT2).AND.(ALT2.GT.ALT3)) GO TO 819
IF((ALT1.GT.ALT2).AND.(ALT2.LT.ALT3)) GO TO 819
IF((ALT4.LT.ALT2).AND.(ALT2.GT.ALT5)) GO TO 819
IF((ALT4.GT.ALT2).AND.(ALT2.LT.ALT5)) GO TO 819
GO TO 99
819 CONTINUE
IAV=IA+IT
IF(IAV.LE.0) GO TO 98
IF(IAV.GT.NX) GO TO 98
JAV=JA+JT
IF(JAV.LE.0) GO TO 98
IF(JAV.GT.NY) GO TO 98
VA1=ALT(IAV,JAV)
UD1=VA1-HV
VD1=VA1-HV
IF(VD1) 1,1,2
C ****
1 IREP=IREP+1
RT=(UA1-HV)/(UA1-VA1)
ALX=HX*(FLOAT(IA-1)
& +RT*FLOAT(IT))
ALY=HY*(FLOAT(JA-1)
& +RT*FLOAT(JT))
CALL TP2TO3(ALX,ALY,2)
ITV=IT
IT=-KS*JT
JT=KS*ITV
UA2=VA1
IF(IREP.GE.4) GO TO 98
GO TO 10
C ****
2 IREP=0
ITV=KS*JT
JTV=-KS*IT
IAV=IAV+ITV
JAV=JAV+JTV
VA2=ALT(IAV,JAV)
IF(VA2-HV) 3,3,4
C ****
3 RT=(VA1-HV)/(VA1-VA2)
ALX=HX*(FLOAT(IA+IT-1)
& +RT*FLOAT(ITV))
ALY=HY*(FLOAT(JA+JT-1)
& +RT*FLOAT(JTV))
CALL TP2TO3(ALX,ALY,2)
KSW(IA,JA)=KV
IA=IA+IT
JA=JA+JT
UA1=VA1
UA2=VA2
GO TO 10
C ****
4 UD2=UA2-HV
ITV=KS*JT
JTV=-KS*IT
RT=(UA2-HV)/(UA2-VA2)
ALX=HX*(FLOAT(IA+ITV-1)
& +RT*FLOAT(IT))
ALY=HY*(FLOAT(JA+JTV-1)
& +RT*FLOAT(JT))
CALL TP2TO3(ALX,ALY,2)
KSW(IA,JA)=KV
IA=IA+IT+ITV
JA=JA+JT+JTV
IT=ITV
JT=JTV
UA1=VA2
C ****
10 CONTINUE
C ****
GO TO 99
98 KSW(IA,JA)=KV
99 CONTINUE
BX=ALX
BY=ALY
RETURN
END
C SUBROUTINE TPCROS ====*====3====*====4====*====5====*====6====*====7
C
C GRID PLOT ROUTINE IN CONT.
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPCROS(NCX,NCY)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/W1(9),IUOD,IHL,IPR,IMOD
COMMON/TPERS3/ALT(2402,2402),HX,HY,NX,NY,HV
COMMON/TPERS6/DFST,DMIN,DMAX,DFAC,DVID
DDFST=DFST
DFST=0.05
XWR=HX*FLOAT(NX-1)
YWR=HY*FLOAT(NY-1)
EEF=HV
IEOH=IUOD
IF(NCX.LE.0) GO TO 30
IYS=NY
DO 10 KX=1,NX,NCX
IF(IYS.EQ.NY) GO TO 12
IYS=NY
YA1=YWR
GO TO 13
12 IYS=1
YA1=0.0
13 XX=HX*FLOAT(KX-1)
IP=3
UA1=ALT(KX,IYS)
IF(UA1.GT.EEF) GO TO 181
IF(IEOH.EQ.-1) GO TO 18
GO TO 182
181 IF(IEOH.EQ.1) GO TO 18
182 CONTINUE
CALL TP2TO3(XX,YA1,IP)
IP=2
CALL TP2TO3(XX,YA1,IP)
18 DO 14 IKY=2,NY
KY=IKY
IF(IYS.EQ.NY) KY=NY-IKY+1
YA2=HY*FLOAT(KY-1)
UA2=ALT(KX,KY)
IF(IP.EQ.3) GO TO 15
IF(UA2.LE.EEF) GO TO 161
IF(IEOH.EQ.-1) GO TO 16
GO TO 162
161 IF(IEOH.EQ.1) GO TO 16
162 CONTINUE
YY=YA1+(YA2-YA1)*((EEF-UA1)/(UA2-UA1))
CALL TP2TO3(XX,YY,IP)
IP=3
GO TO 17
16 CALL TP2TO3(XX,YA2,IP)
GO TO 17
15 IF(UA2.GT.EEF) GO TO 171
IF(IEOH.EQ.-1) GO TO 17
GO TO 172
171 IF(IEOH.EQ.1) GO TO 17
172 CONTINUE
YY=YA1+(YA2-YA1)*((EEF-UA1)/(UA2-UA1))
CALL TP2TO3(XX,YY,IP)
IP=2
CALL TP2TO3(XX,YY,IP)
CALL TP2TO3(XX,YA2,IP)
17 UA1=UA2
YA1=YA2
14 CONTINUE
10 CONTINUE
30 IF(NCY.LE.0) GO TO 31
IXS=1
DO 20 KY=1,NY,NCY
IF(IXS.EQ.NX) GO TO 22
IXS=NX
XA1=XWR
GO TO 23
22 IXS=1
XA1=0.0
23 YY=HY*FLOAT(KY-1)
IP=3
UA1=ALT(IXS,KY)
IF(UA1.GT.EEF) GO TO 281
IF(IEOH.EQ.-1) GO TO 28
GO TO 282
281 IF(IEOH.EQ.1) GO TO 28
282 CONTINUE
CALL TP2TO3(XA1,YY,IP)
IP=2
CALL TP2TO3(XA1,YY,IP)
28 DO 24 IKX=2,NX
KX=IKX
IF(IXS.EQ.NX) KX=NX-IKX+1
XA2=HX*FLOAT(KX-1)
UA2=ALT(KX,KY)
IF(IP.EQ.3) GO TO 25
IF(UA2.LE.EEF) GO TO 261
IF(IEOH.EQ.-1) GO TO 26
GO TO 262
261 IF(IEOH.EQ.1) GO TO 26
262 CONTINUE
XX=XA1+(XA2-XA1)*((EEF-UA1)/(UA2-UA1))
CALL TP2TO3(XX,YY,IP)
IP=3
GO TO 27
26 CALL TP2TO3(XA2,YY,IP)
GO TO 27
25 IF(UA2.GT.EEF) GO TO 271
IF(IEOH.EQ.-1) GO TO 27
GO TO 272
271 IF(IEOH.EQ.1) GO TO 27
272 CONTINUE
XX=XA1+(XA2-XA1)*((EEF-UA1)/(UA2-UA1))
CALL TP2TO3(XX,YY,IP)
IP=2
CALL TP2TO3(XX,YY,IP)
CALL TP2TO3(XA2,YY,IP)
27 UA1=UA2
XA1=XA2
24 CONTINUE
20 CONTINUE
31 DFST=DDFST
RETURN
END
C SUBROUTINE TP2TO3 ====*====3====*====4====*====5====*====6====*====7
C
C WHEN IHL.GE.0 THEN
C TRANSFORMATION FROM TWO DIMENSIONAL PLANE TO
C THEREE DIMENSIONAL COORDINATES
C WHEN IHL.LT.0 THEN
C DIRECTLY CALL PLOTPS
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TP2TO3(ALX,ALY,IP)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 KB
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS4/XOR(3),WW(3,3),XFC,YFC
DIMENSION XD(3),KB(3)
C WRITE(6,*) ALX,ALY,XOR
JP=1
IF(IHL.EQ.-2) GO TO 811
XD(1)=ALX
XD(2)=ALY
DO 71 I1=1,3
WA=XOR(I1)
DO 72 I2=1,2
72 WA=WA+WW(I1,I2)*XD(I2)
IF(IMOD.EQ.1) THEN
IF(WA.GT.XM(I1,2)) JP=0
IF(WA.LT.XM(I1,1)) JP=0
ELSE IF(IMOD.EQ.-1) THEN
IF(WA.GT.XM(I1,3)) JP=0
IF(WA.LT.0.0) JP=0
END IF
71 KB(I1)=WA
C WRITE(6,*) KB
IF(IHL.LT.0) GO TO 811
CALL TPPLOT(KB(1),KB(2),KB(3),IP)
RETURN
811 IF(JP.EQ.0) RETURN
IF(IMOD.EQ.1) CALL PLOTPS(ALX*XFC,ALY*YFC,IP-1)
IF(IMOD.EQ.-1) CALL PLOTPS(ALX,ALY,IP-1)
RETURN
END
SUBROUTINE PLOTPS(AX,AY,IP)
IMPLICIT REAL*8(A-H,O-Z)
REAL*4 XX,YY
XX=AX
YY=AY
IF(IP.EQ.2) THEN
CALL MOVETO(XX,YY)
ELSE IF(IP.EQ.1) THEN
CALL LINETO(XX,YY)
ELSE
WRITE(6,*) ' STOP IN PLOTPS IP=',IP
END IF
RETURN
END
C SUBROUTINE TPLIN2 ====*====3====*====4====*====5====*====6====*====7
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPLIN2(X,Y,IP)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS0/GENT(2)
COMMON/TPERS9/XOM(3),WM(3,2)
DIMENSION XD(2),XA(3)
JP=IP+1
XD(1)=GENT(1)+X
XD(2)=GENT(2)+Y
DO 771 I1=1,3
WA=XOM(I1)
DO 772 I2=1,2
772 WA=WA+WM(I1,I2)*XD(I2)
771 XA(I1)=WA
CALL TPPLOT(XA(1),XA(2),XA(3),JP)
RETURN
END
SUBROUTINE TPORIG(XO,YO)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS0/GENT(2)
GENT(1)=GENT(1)+XO
GENT(2)=GENT(2)+YO
RETURN
END
CC SUBROUTINE TPCIRC ====*====3====*====4====*====5====*====6====*====7
C
C CIRCLE
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPCIRC(XC,YC,RR)
IMPLICIT REAL*8(A-H,O-Z)
C write(6,*) ' tpcirc',xc,yc,rr
CALL TPORIG(XC,YC)
CALL TPLIN2(RR,0.0D0,2)
DO 1 I=1,72
TH=FLOAT(I*5)*0.01745329
X=RR*COS(TH)
Y=RR*SIN(TH)
CALL TPLIN2(X,Y,1)
1 CONTINUE
CALL TPORIG(-XC,-YC)
RETURN
END
C SUBROUTINE TPLANE ====*====3====*====4====*====5====*====6====*====7
C
C PLANE IS SET WITH X-AXIS IN DX DIRECTION
C Y-AXIS IN DY DIRECTION
C AND ORIGIN AT C
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPLANE(DX,DY,C)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS0/GENT(2)
COMMON/TPERS9/XOM(3),WM(3,2)
DIMENSION DX(3),DY(3),C(3),AX(3),AY(3)
WX=SQRT(DX(1)*DX(1)+DX(2)*DX(2)+DX(3)*DX(3))
IF(WX.GT.1.0D-7) GO TO 1
WRITE(6,*) ' STOP IN TPLANE, NORM OF DX=',WX
STOP
1 AX(1)=DX(1)/WX
AX(2)=DX(2)/WX
AX(3)=DX(3)/WX
WY=SQRT(DY(1)*DY(1)+DY(2)*DY(2)+DY(3)*DY(3))
IF(WY.GT.10D-7) GO TO 2
WRITE(6,*) ' STOP IN TPLANE, NORM OF DY=',WX
STOP
2 AY(1)=DY(1)/WY
AY(2)=DY(2)/WY
AY(3)=DY(3)/WY
DO 4 I=1,3
WM(I,1)=AX(I)
WM(I,2)=AY(I)
XOM(I)=C(I)
4 CONTINUE
DO 5 I=1,2
GENT(I)=0.0
5 CONTINUE
RETURN
END
C SUBROUTINE TPCLEA ====*====3====*====4====*====5====*====6====*====7
C
C MAIN ROUTINE FOR PLOT OF SPHERES AND BARS
C
C 1989/12/09 BY A.YANASE
C 2002/09/26 BY A.YANASE CHANG EENTRY TO SUBROUTINE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPCLEA
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),XB(10,300),RB(300),
&IB(2,300),NS,NB
C EXTERNAL TPMASB
CALL TPHILP(0.0D0)
CALL TPCIHL(1)
NS=0
NB=0
RETURN
END
C SUBROUTINE TPSETS ====*====3====*====4====*====5====*====6====*====7
C
C SET SPHERE
C 2002/09/26 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPSETS(XX,R,IP)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),XB(10,300),RB(300),
&IB(2,300),NS,NB
DIMENSION XX(3)
IF(IP.GT.9.OR.IP.LE.-1) RETURN
NS=NS+1
IF(IMOD.EQ.-1) GO TO 2
DO 1 I=1,3
1 XC(I,NS)=XX(I)
RR(NS)=R
IIP(NS)=IP
RETURN
2 DO 3 I=1,3
3 XC(I,NS)=XX(I)/FCTR(I)+XM(I,1)
RR(NS)=R/FCTR(1)
IIP(NS)=IP
RETURN
END
C SUBROUTINE TPSETB ====*====3====*====4====*====5====*====6====*====7
C
C SET BAR
C 2002/09/26 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPSETB(N1,N2,RL)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),XB(10,300),RB(300),
&IB(2,300),NS,NB
DIMENSION DX(3),DY(3),EY(3)
IF(N1.GT.NS) RETURN
IF(N2.GT.NS) RETURN
NB=NB+1
WA=0.0
DO 11 I=1,3
XB(I,NB)=XC(I,N1)-XC(I,N2)
11 WA=WA+XB(I,NB)**2
XB(4,NB)=WA
IF(IMOD.EQ.1) RB(NB)=RL
IF(IMOD.EQ.-1) RB(NB)=(RL/FCTR(1))
DO 12 I=1,3
12 EY(I)=(EYE(I)/FCTR(I))+XM(I,1)
DX(1)=XB(2,NB)*(EY(3)-XC(3,N1))-XB(3,NB)*(EY(2)-XC(2,N1))
DX(2)=XB(3,NB)*(EY(1)-XC(1,N1))-XB(1,NB)*(EY(3)-XC(3,N1))
DX(3)=XB(1,NB)*(EY(2)-XC(2,N1))-XB(2,NB)*(EY(1)-XC(1,N1))
DY(1)=XB(2,NB)*DX(3)-XB(3,NB)*DX(2)
DY(2)=XB(3,NB)*DX(1)-XB(1,NB)*DX(3)
DY(3)=XB(1,NB)*DX(2)-XB(2,NB)*DX(1)
WC=0.0
WD=0.0
DO 13 I=1,3
WC=WC+DX(I)*DX(I)
13 WD=WD+DY(I)*DY(I)
WC=SQRT(WC)
IF(WC.LT.1.E-7) GO TO 16
WD=SQRT(WD)
DO 18 I=1,3
DX(I)=DX(I)/WC
18 DY(I)=DY(I)/WD
WE=0.0
DO 14 I=1,2
DO 15 J=I+1,3
WE=WE+((EY(J)-XC(J,N1))*XB(I,NB)- (EY(I)-XC(I,N1))*XB(J,NB))**2
15 CONTINUE
14 CONTINUE
A=SQRT(WE/XB(4,NB))
P=(RL*RL)/A
RM=RL*SQRT(1.0-P/A)
DO 17 I=1,3
XB(4+I,NB)=P*DY(I)+RM*DX(I)
XB(7+I,NB)=P*DY(I)-RM*DX(I)
17 CONTINUE
IB(1,NB)=N1
IB(2,NB)=N2
RETURN
16 NB=NB-1
RETURN
END
C SUBROUTINE TPDRWS ====*====3====*====4====*====5====*====6====*====7
C
C DRAW SPHERE
C 2002/09/26 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPDRWS(N)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),XB(10,300),RB(300),
&IB(2,300),NS,NB
IF(N.GT.NS) RETURN
CALL PTIME(TB)
IM=IMOD
IMOD=1
X=XC(1,N)
Y=XC(2,N)
Z=XC(3,N)
RS=RR(N)
C write(6,*) ' tpdrws',n,rs,iip(n)
IF(IIP(N).EQ.0) CALL TPSPHE(X,Y,Z,RS,0,0)
IF(IIP(N).EQ.1) CALL TPSPHR(X,Y,Z,RS,1)
IF(IIP(N).EQ.2) CALL TPSPHE(X,Y,Z,RS,4,3)
IF(IIP(N).EQ.3) CALL TPSPHR(X,Y,Z,RS,3)
IF(IIP(N).EQ.4) CALL TPSPHE(X,Y,Z,RS,8,7)
IF(IIP(N).EQ.5) CALL TPSPHR(X,Y,Z,RS,5)
IF(IIP(N).EQ.6) CALL TPSPHE(X,Y,Z,RS,10,9)
IF(IIP(N).EQ.7) CALL TPSPHR(X,Y,Z,RS,7)
IF(IIP(N).EQ.8) CALL TPSPHE(X,Y,Z,RS,12,11)
IF(IIP(N).EQ.9) CALL TPSPHR(X,Y,Z,RS,9)
IMOD=IM
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1.AND.MOD(IIP(N),2).EQ.0) WRITE(6,6003) TE,X,Y,Z,RS
6003 FORMAT(' TIME FOR TPSPHE TPDRWS',F7.3,4E11.4)
IF(IPR.GE.1.AND.MOD(IIP(N),2).EQ.1) WRITE(6,6007) TE,X,Y,Z,RS
6007 FORMAT(' TIME FOR TPSPHR TPDRWS',F7.3,4E11.4)
RETURN
END
C SUBROUTINE TPDRWB ====*====3====*====4====*====5====*====6====*====7
C
C DRAW BAR
C 2002/09/26 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPDRWB(M)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),XB(10,300),RB(300),
&IB(2,300),NS,NB
DIMENSION XA(3),XAA(3)
IF(M.GT.NB) RETURN
CALL PTIME(TB)
IM=IMOD
IMOD=1
N1=IB(1,M)
N2=IB(2,M)
RQ=RB(M)
CALL TPINSS(N1,N2,RQ)
DO 22 J=1,2
DO 21 I=1,3
II=I+(J-1)*3+4
XA(I)=XC(I,N1)+XB(II,M)
XAA(I)=XC(I,N2)+XB(II,M)
21 CONTINUE
CALL TPLINE(XA(1),XA(2),XA(3),XAA(1),XAA(2),XAA(3),30)
22 CONTINUE
IMOD=IM
CALL PTIME(TA)
TE=(TA-TB)*60.0*60.0
IF(IPR.GE.1) WRITE(6,6001) TE,(XC(I,N1),I=1,3)
& ,(XC(I,N2),I=1,3)
6001 FORMAT(' TIME TPLINE TPDRWB',F7.4,3E11.4/26X,3E11.4)
RETURN
END
C SUBROUTINE TPLINS ====*====3====*====4====*====5====*====6====*====7
C
C DRAW LINE BETWEEN SHERES
C 2002/09/26 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPLINS(N1,N2)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),XB(10,300),RB(300),
&IB(2,300),NS,NB
IF(N1.GT.NS) RETURN
IF(N2.GT.NS) RETURN
X1=XC(1,N1)
Y1=XC(2,N1)
Z1=XC(3,N1)
X2=XC(1,N2)
Y2=XC(2,N2)
Z2=XC(3,N2)
CALL TPLINE(X1,Y1,Z1,X2,Y2,Z2,60)
RETURN
END
C SUBROUTINE TPSPHE ====*====3====*====4====*====5====*====6====*====7
C
C PLOT ROUTINE OF SPHERE WITH PATER OF EARTH GLOVE
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPSPHE(XC,YC,ZC,RS,NPHI,NTHITA)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
DIMENSION DX(3),DY(3),C(3),X(3),PP(3)
C write(6,*) ' tpsphe',xc,yc,zc,rs,nphi,nthita
IM=IMOD
IF(IMOD.EQ.-1) GO TO 1
IMOD=-1
X(1)=(XC-XM(1,1))*FCTR(1)
X(2)=(YC-XM(2,1))*FCTR(2)
X(3)=(ZC-XM(3,1))*FCTR(3)
RR=RS*FCTR(1)
GO TO 2
1 X(1)=XC
X(2)=YC
X(3)=ZC
RR=RS
2 NP=NPHI
NT=NTHITA
WW=0.0
DO 3 I=1,3
PP(I)=EYE(I)-X(I)
WW=WW+PP(I)*PP(I)
3 CONTINUE
A=SQRT(WW)
IF(A-RR.LT.0.01) GO TO 7
PX=PP(1)/A
PY=PP(2)/A
PZ=PP(3)/A
P=(RR*RR)/A
R=RR*SQRT(1.0-P/A)
WA=SQRT(PX*PX+PY*PY)
DX(1)=PY/WA
DX(2)=-PX/WA
DX(3)=0.0
DY(1)=(PZ*PX)/WA
DY(2)=(PZ*PY)/WA
DY(3)=-WA
C(1)=X(1)+P*PX
C(2)=X(2)+P*PY
C(3)=X(3)+P*PZ
CALL TPLANE(DX,DY,C)
CALL TPCIRC(0.0D0,0.0D0,R)
IF(NP.LE.0) GO TO 4
DO 5 I=1,3
C(I)=X(I)
5 DX(I)=0.0
DX(3)=1.0
DY(3)=0.0
W1=3.141592/FLOAT(NP)
DO 6 I=1,NP
PHAI=FLOAT(I-1)*W1
DY(1)=COS(PHAI)
DY(2)=SIN(PHAI)
CALL TPLANE(DX,DY,C)
CALL TPCIRC(0.0D0,0.0D0,RR)
6 CONTINUE
4 IF(NT.LE.0) GO TO 7
DO 8 I=1,3
DX(I)=0.0
DY(I)=0.0
8 C(I)=X(I)
DX(1)=1.0
DY(2)=1.0
W2=3.141592/FLOAT(NT+1)
DO 9 I=1,NT
THI=FLOAT(I)*W2
C(3)=X(3)+RR*COS(THI)
CALL TPLANE(DX,DY,C)
RA=RR*SIN(THI)
CALL TPCIRC(0.0D0,0.0D0,RA)
9 CONTINUE
7 DO 10 I=1,3
DX(I)=0.0
DY(I)=0.0
10 C(I)=0.0
DX(1)=1.0
DY(2)=1.0
CALL TPLANE(DX,DY,C)
IF(IM.EQ.-1) GO TO 11
CALL TPUNIT(IM)
11 RETURN
END
C SUBROUTINE TPSPHE ====*====3====*====4====*====5====*====6====*====7
C
C PLOT ROUTINE OF SPHERE WITH PATERN OF CIRCLE ON PLANES
C PERPENDICLAR TO THREE DIRECTION OF X,Y AND Z
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPSPHR(XC,YC,ZC,RS,NZ)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(8),EYE(3)
DIMENSION DX(3),DY(3),C(3),X(3),PP(3)
IM=IMOD
IF(IMOD.EQ.-1) GO TO 1
IMOD=-1
X(1)=(XC-XM(1,1))*FCTR(1)
X(2)=(YC-XM(2,1))*FCTR(2)
X(3)=(ZC-XM(3,1))*FCTR(3)
RR=RS*FCTR(1)
GO TO 2
1 X(1)=XC
X(2)=YC
X(3)=ZC
RR=RS
2 N=NZ
WW=0.0
DO 3 I=1,3
PP(I)=EYE(I)-X(I)
WW=WW+PP(I)*PP(I)
3 CONTINUE
A=SQRT(WW)
IF(A-RR.LT.0.01) GO TO 7
PX=PP(1)/A
PY=PP(2)/A
PZ=PP(3)/A
P=(RR*RR)/A
R=RR*SQRT(1.0-P/A)
WA=SQRT(PX*PX+PY*PY)
DX(1)=PY/WA
DX(2)=-PX/WA
DX(3)=0.0
DY(1)=(PZ*PX)/WA
DY(2)=(PZ*PY)/WA
DY(3)=-WA
C(1)=X(1)+P*PX
C(2)=X(2)+P*PY
C(3)=X(3)+P*PZ
CALL TPLANE(DX,DY,C)
CALL TPCIRC(0.0D0,0.0D0,R)
IF(N.LE.0) GO TO 7
DO 21 II=1,2
DO 22 JJ=II+1,3
KK=JJ-II
IF(JJ.EQ.2) KK=3
DO 24 I=1,3
C(I)=X(I)
DX(I)=0.0
24 DY(I)=0.0
DX(II)=1.0
DY(JJ)=1.0
DZ=(2.0*RR)/FLOAT(N+1)
DO 23 I=1,N
RB=-RR+DZ*FLOAT(I)
RA=SQRT(RR*RR-RB*RB)
C(KK)=RB+X(KK)
CALL TPLANE(DX,DY,C)
CALL TPCIRC(0.0D0,0.0D0,RA)
23 CONTINUE
22 CONTINUE
21 CONTINUE
7 DO 10 I=1,3
DX(I)=0.0
DY(I)=0.0
10 C(I)=0.0
DX(1)=1.0
DY(2)=1.0
CALL TPLANE(DX,DY,C)
IF(IM.EQ.-1) GO TO 11
CALL TPUNIT(IM)
11 RETURN
END
C SUBROUTINE TPINSS ====*====3====*====4====*====5====*====6====*====7
C
C PLOT ROUTINE OF CROSS SECTION SPHERE AND BAR
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE TPINSS(N1,N2,R)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/XM(3,3),IUOD,IHL,IPR,IMOD
COMMON/TPERS2/FCTR(3),W(11)
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),
& XB(10,300),RB(300),IB(2,300),NS,NB
DIMENSION X1(3),X2(3),DELT(3),CENT1(3),CENT2(3), X(3),Y(3),DX(3)
IM=IMOD
IMOD=-1
RF=R*FCTR(1)
R1=RR(N1)*FCTR(1)
R2=RR(N2)*FCTR(1)
DO 10 I=1,3
X1(I)=(XC(I,N1)-XM(I,1))*FCTR(I)
X2(I)=(XC(I,N2)-XM(I,1))*FCTR(I)
DELT(I)=X2(I)-X1(I)
10 CONTINUE
IF(R1.LE.RF) GO TO 11
RC1=SQRT(R1*R1-RF*RF)
GO TO 12
11 RC1=0.0
12 CONTINUE
IF(R2.LE.RF) GO TO 13
RC2=SQRT(R2*R2-RF*RF)
GO TO 14
13 RC2=0.0
14 CONTINUE
ANORM=SQRT(DELT(1)*DELT(1)+DELT(2)*DELT(2)+DELT(3)*DELT(3))
DO 20 I=1,3
DX(I)=DELT(I)/ANORM
CENT1(I)=X1(I)+RC1*DX(I)
CENT2(I)=X2(I)-RC2*DX(I)
20 CONTINUE
IF(ABS(1.0-ABS(DX(3))).GT.1.0E-4) GO TO 30
X(1)=1.0
X(2)=0.0
X(3)=0.0
Y(1)=0.0
Y(2)=1.0
Y(3)=0.0
GO TO 31
30 CONTINUE
X(1)=DX(2)
X(2)=-DX(1)
X(3)=0.0
Y(1)=-DX(1)*DX(3)
Y(2)=-DX(2)*DX(3)
Y(3)=DX(2)*DX(2)+DX(1)*DX(1)
31 CONTINUE
CALL TPLANE(X,Y,CENT1)
CALL TPCIRC(0.0D0,0.0D0,RF)
CALL TPLANE(X,Y,CENT2)
CALL TPCIRC(0.0D0,0.0D0,RF)
IMOD=IM
RETURN
END
C FUNCTION TPMASB 2====*====3====*====4====*====5====*====6====*====7
C
C HIDDEN LINE FUCTION FOR SPHERES AND BARS PLOT
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
FUNCTION TPMASB(Q,P)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TPERS1/X0,Y0,Z0,XN,YN,ZN,XW,YW,ZW,IUOD,IHL,IPR,IMOD
COMMON/TPERS2/XFCTR,YFCTR,ZFCTR,XP0,CA,SA,ET,EL,SG,ZM,ZC
& ,XEYE,YEYE,ZEYE
COMMON/TPERS7/XC(3,300),RR(300),IIP(300),XB(10,300), RB(300),
&IB(2,300),NS,NB
DIMENSION Q(3),P(3),X(3),X1(3),X2(3),DELT(3),EYE(3)
EYE(1)=XEYE/XFCTR+X0
EYE(2)=YEYE/YFCTR+Y0
EYE(3)=ZEYE/ZFCTR+Z0
IF(NS.EQ.0) GO TO 11
DO 10 I=1,NS
DO 3 J=1,3
3 X(J)=XC(J,I)
R=RR(I)
DO 5 IJ=1,3
IF(Q(IJ).LT.X(IJ)-R.AND.EYE(IJ).LT.X(IJ)-R) GO TO 10
IF(Q(IJ).GT.X(IJ)+R.AND.EYE(IJ).GT.X(IJ)+R) GO TO 10
5 CONTINUE
R=R*0.999
SL=TPSLEN(Q,P,X,R)
IF(SL.LT.0.0) GO TO 100
10 CONTINUE
11 CONTINUE
IF(NB.EQ.0) GO TO 21
DO 20 I=1,NB
IB1=IB(1,I)
IB2=IB(2,I)
DO 13 J=1,3
X1(J)=XC(J,IB1)
13 X2(J)=XC(J,IB2)
DO 15 IJ=1,3
XMIN=DMIN1(X1(IJ),X2(IJ))-RB(I)
XMAX=DMAX1(X1(IJ),X2(IJ))+RB(I)
IF(Q(IJ).LT.XMIN.AND.EYE(IJ).LT.XMIN) GO TO 20
IF(Q(IJ).GT.XMAX.AND.EYE(IJ).GT.XMAX) GO TO 20
15 CONTINUE
R=RB(I)
R=R*0.999
DELT(1)=XB(1,I)
DELT(2)=XB(2,I)
DELT(3)=XB(3,I)
BL=TPBLEN(Q,P,X1,X2,DELT,XB(4,I),R)
IF(BL.LT.0.0) GO TO 100
20 CONTINUE
21 CONTINUE
TPMASB=1.0
RETURN
100 TPMASB=-1.0
RETURN
END
C FUNCTION TPSLEN 2====*====3====*====4====*====5====*====6====*====7
C
C HIDDEN LINE FUCTION FOR SPHERES
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
FUNCTION TPSLEN(Q,P,X,R)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION Q(3),P(3),X(3)
DELT2=0.0
DO 10 I=1,2
DO 10 J=I+1,3
DELT2=DELT2+((X(J)-Q(J))*P(I)-(X(I)-Q(I))*P(J))**2
10 CONTINUE
TPSLEN=SQRT(DELT2)-R
IF(TPSLEN.LT.0.0) GO TO 1
TPSLEN=1.0
RETURN
1 CONTINUE
A=0.0
DO 3 I=1,3
A=A+(Q(I)-X(I))*(Q(I)-X(I))
3 CONTINUE
IF(A.LT.R*R) GO TO 99
B=0.0
DO 13 I=1,3
B=B+(Q(I)-X(I))*P(I)
13 CONTINUE
IF(B.LT.0.0) GO TO 99
TPSLEN=1.0
RETURN
99 TPSLEN=-1.0
RETURN
END
C FUNCTION TPBLEN 2====*====3====*====4====*====5====*====6====*====7
C
C HIDDEN LINE FUCTION FOR BARS
C
C 1989/12/09 BY A.YANASE
C
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
FUNCTION TPBLEN(Q,P,X,Y,DELT,ANOR,R)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION Q(3),P(3),X(3),Y(3),DELT(3)
Q1=Q(1)
Q2=Q(2)
Q3=Q(3)
P1=P(1)
P2=P(2)
P3=P(3)
X1=X(1)
X2=X(2)
X3=X(3)
Y1=Y(1)
Y2=Y(2)
Y3=Y(3)
D1=DELT(1)
D2=DELT(2)
D3=DELT(3)
A1=D2*P3-P2*D3
A2=D3*P1-P3*D1
A3=D1*P2-P1*D2
XQ1=X1-Q1
XQ2=X2-Q2
XQ3=X3-Q3
A=ABS(XQ1*A1+XQ2*A2+XQ3*A3)
B=SQRT(A1*A1+A2*A2+A3*A3)
TP=A/B-R
IF(TP.GE.0.0) GO TO 6
AS=A1*(XQ3*P2-XQ2*P3)
& +A2*(XQ1*P3-XQ3*P1)
& +A3*(XQ2*P1-XQ1*P2)
YQ1=Y1-Q1
YQ2=Y2-Q2
YQ3=Y3-Q3
AE=A1*(YQ3*P2-YQ2*P3)
& +A2*(YQ1*P3-YQ3*P1)
& +A3*(YQ2*P1-YQ1*P2)
IF(AS*AE.GE.0.0) GO TO 2
BS=A1*(XQ2*D3-XQ3*D2)
& +A2*(XQ3*D1-XQ1*D3)
& +A3*(XQ1*D2-XQ2*D1)
IF(BS.LT.0.0) GO TO 7
TB=(XQ1*D2-XQ2*D1)**2
& +(XQ2*D3-XQ3*D2)**2
& +(XQ3*D1-XQ1*D3)**2
TB=SQRT(TB/ANOR)-R
IF(TB.LT.0.0) GO TO 7
IF(TB.GE.0.0) GO TO 6
2 PD=P1*D1+P2*D2+P3*D3
IF(ABS(PD).LE.1.0E-5) GO TO 6
IF(AE.GE.0.0) GO TO 1
XQ1=YQ1
XQ2=YQ2
XQ3=YQ3
1 XQD=D1*XQ1+D2*XQ2+D3*XQ3
IF(XQD*PD.LT.1.0E-5) GO TO 6
T=XQD/PD
TL=SQRT((XQ1-P1*T)**2
& +(XQ2-P2*T)**2
& +(XQ3-P3*T)**2)
IF(TL.GE.R) GO TO 6
7 TPBLEN=-1.0
RETURN
6 TPBLEN=1.0
RETURN
END