home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e003 / 12.ddi / CONY.FOR next >
Encoding:
Text File  |  1988-01-04  |  13.4 KB  |  471 lines
  1. C     SUBROUTINE CONY  CONECT SAP6 WITH TAP6  IT IS USED TO GENERATION  
  2. C     TEMPRATURE TO SAP6 AND THROUGH SAP6 CALCULAT THERMAL STREE AND
  3. C     DISPLACEMENT. CONY CAN READ SAP6 DATA CARD AND TRANFORM TO TAP6 DIF
  4. C     ERECE NETWORK. WE GET THIS VERSION SUB CONY HAS MANY ERROR SO CAN'T
  5. C     BE USED.
  6. C      THIS SUB IS MODIFIED BY T.U.P  CAD OFFICE IN BEIJING, 1986
  7. C
  8.       SUBROUTINE CONY
  9.       INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000),
  10.      X        NXP(1000),KEY(1000),IZ(1000),IPR(500),IPH(100),CON(20),
  11.      X        NTHEAD(100),NTTAIL(100),NFLAG(1000),MBEAM(100),MTRUS(100)
  12.       DIMENSION ZZ(16000),DYTCQ(11030),Y(3000),T(1000),C(1000),Q(1000),
  13.      X          GAIN(500),A(1000),B(1000),D(1000),RHO(50),CP(50),TK(50),
  14.      X          TPROG(2),AA(7),IX(8),LX(8),XX(3),XO(3),KS(4),THK(100),
  15.      X          BSECA(100),TSECA(100),AREA(3),ANGL(3)
  16.       REAL*8 TEST,TEST1,TEST2,TEST3,TEST4,TEST5,TEST6,TEST7,TEST8,TEST9
  17.       COMMON /I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IRP,IPH,NTHEAD,
  18.      X        NTTAIL,NFLAG,IZ
  19.       COMMON /LAB1/PLOT(100,10),IPLOT(10)
  20.       COMMON /LAB2/DATE(2)
  21.       COMMON /LAB3/ICASE,XLIM(2),YLIM(2),CT(5),XLAB(5),YLAB(5)
  22.       COMMON /LAB4/SPROG(2)
  23.       COMMON /LAB5/TITLE
  24.       COMMON ZZ
  25.       EQUIVALENCE (ZZ(1),TIME),(ZZ(2),Y(1),DYTCQ(1)),
  26.      X    (ZZ(3002),T(1)),(ZZ(4002),C(1)),(ZZ(11061),MAX1T)
  27.      X           ,(ZZ(5002),Q(1)),(ZZ(6002),D(1)),(ZZ(11060),MIN1T),
  28.      X    (ZZ(7002),A(1)),(ZZ(8002),B(1)),(ZZ(11077),IMIN),
  29.      X    (ZZ(11078),IMAX),(ZZ(11085),MAXNOT),(ZZ(11086),MAXNOY),
  30.      X    (ZZ(11002),CON(1))
  31.       DATA TEST1/8HELEMENTS/
  32.       DATA TEST2/8HNODE-INP/
  33.       DATA TEST3/8HMATERIAL/
  34.       DATA TEST4/8HEXECUTE-/
  35.       DATA TEST5/8HREADDATA/
  36.       DATA TEST6/8HTHIK-TYP/
  37.       DATA TEST7/8HTRUS-TYP/
  38.       DATA TEST8/8HBEAM-TYP/
  39.       DATA TEST9/8HTEMPDIST/
  40.       DATA KS/1H ,1HR,1HC,1HS/
  41.       REWIND 19
  42.       M=1000
  43.       PI2=(3.14159/180.0)**2
  44. C
  45. C     SEAARCH THE SECTION HEAD WORD 'NODE-INP'
  46. C
  47. 10    READ(19,900) TEST
  48.       IF(TEST.NE.TEST2) GO TO 10
  49.       DO 2  J=1,8
  50.       IX(J)=0
  51. 2     CONTINUE
  52.           DO 134 JJ=1,3
  53. 134       XX(JJ)=0.0
  54. 11    DO 5 J=1,3
  55.       XO(J)=XX(J)
  56. 5     LX(J)=IX(J)
  57. C
  58. C     BEGIN TO READ NODE INFORMATION
  59. C
  60.       READ(19,920) IX(1),XX(1),XX(2),XX(3),IX(2),KT,IX(3),XY
  61.       IF(IX(1).LE.0) GO TO 20
  62.       L2=IX(1)
  63.       MAX1T=MAX0(MAX1T,L2)
  64.       MIN1T=MIN0(MIN1T,L2)
  65.       MAXNOT=MAX1T
  66.       IF(KT.EQ.KS(3)) GO TO 6
  67.       IF(KT.EQ.KS(4)) GO TO 7
  68.       NFLAG(L2)=2
  69.       A(L2)=XX(1)
  70.       B(L2)=XX(2)
  71.       D(L2)=XX(3)
  72.       GO TO 8
  73. 6     NFLAG(L2)=3
  74.       PHI=XX(2)
  75.       A(L2)=XX(1)*COS(PHI)
  76.       B(L2)=XX(1)*SIN(PHI)
  77.       D(L2)=XX(3)
  78.       GO TO 8
  79. 7     NFLAG(L2)=4
  80.       PHI=XX(2)
  81.       THETA=XX(3)
  82.       A(L2)=XX(1)*SIN(THETA)*COS(PHI)
  83.       B(L2)=XX(1)*SIN(THETA)*SIN(PHI)
  84.       D(L2)=XX(1)*COS(THETA)
  85. 8     IF(IX(3).LE.0) GO TO 11
  86.       L1=LX(1)
  87.       L3=IX(3)
  88.       DN=0.0
  89.       X=0.01*XY
  90.       IF(XY.LE.0.0) X=1.0
  91.       LDX=L2-L1
  92.       K=0
  93.       DO 12 J=1,LDX,L3
  94.       K=K+1
  95. 12    DN=DN+X**(K-1)
  96.       DX1=(XX(1)-XO(1))/DN
  97.       DX2=(XX(2)-XO(2))/DN
  98.       DX3=(XX(3)-XO(3))/DN
  99.       K=0
  100.       DO 15 J=L1,L2,L3
  101.       K=K+1
  102.       IF(J.EQ.1) GO TO 15
  103.       XX(1)=XO(1)+DX1*(K-1)*X**(J-L1)
  104.       XX(2)=XO(2)+DX2*(K-1)*X**(J-L1)
  105.       XX(3)=XO(3)+DX3*(K-1)*X**(J-L1)
  106. 132   IF(KT.EQ.KS(3)) GO TO 16
  107.       IF(KT.EQ.KS(4)) GO TO 17
  108.       NFLAG(J)=2
  109.       A(J)=XX(1)
  110.       B(J)=XX(2)
  111.       D(J)=XX(3)
  112.       GO TO 15
  113. 16    NFLAG(J)=3
  114.       PHI=XX(2)
  115.       A(J)=XX(1)*COS(PHI)
  116.       B(J)=XX(1)*SIN(PHI)
  117.       D(J)=XX(3)
  118.       GO TO 15
  119. 17    NFLAG(J)=4
  120.       PHI=XX(2)
  121.       THETA=XX(3)
  122.       A(J)=XX(1)*SIN(THETA)*COS(PHI)
  123.       B(J)=XX(1)*SIN(THETA)*SIN(PHI)
  124.       D(J)=XX(1)*COS(THETA)
  125. C8    XO(1)=XX(1)
  126. C     XO(2)=XX(2)
  127. C     XO(3)=XX(3)
  128. 15    CONTINUE
  129.       GO TO 11
  130. 20    REWIND 19
  131. C
  132. C     SEARCH THE SECTION HEAD WORD 'TRUS-TYP'
  133. C
  134. 21    READ(19,900) TEST
  135.       IF(TEST.EQ.TEST4.OR.TEST.EQ.TEST5) GO TO 30
  136.       IF(TEST.NE.TEST7) GO TO 21
  137. C
  138. C     BEGIN TO READ TRUS-TYP PROERTIES
  139. C
  140. 22    READ(19,970) M1,M2,XY
  141.       IF(M1.LE.0) GO TO 30
  142.       TSECA(M1)=XY
  143.       MTRUS(M1)=M2
  144.       GO TO 22
  145. 30    REWIND 19
  146. C
  147. C     SEARCH THE SECTION HEAD WORD 'BEAM-TYP'
  148. C
  149. 31    READ(19,900) TEST
  150.       IF(TEST.EQ.TEST4.OR.TEST.EQ.TEST5) GO TO 40
  151.       IF(TEST.NE.TEST8) GO TO 31
  152. C
  153. C     BEGIN TO READ BEAM-TYP PROPERTIES
  154. C
  155. 32    READ(19,970) M1,M2,XY
  156.       IF(M1.LE.0) GO TO 40
  157.       BSECA(M1)=XY
  158.       MBEAM(M1)=M2
  159.       GO TO 32
  160. 40    REWIND 19
  161. C
  162. C     SEARCH THE SECTION HEAD WORD 'THIK-TYP'
  163. C
  164. 41    READ(19,900) TEST
  165.       IF(TEST.EQ.TEST4.OR.TEST.EQ.TEST5) GO TO 50
  166.       IF(TEST.NE.TEST6) GO TO 41
  167. C
  168. C     BEGIN TO READ THE THICKNESS OF PLATE-SHELL EMEMENTS
  169. C
  170. 42    READ(19,930) M1,XY
  171.       IF(M1.LE.0) GO TO 50
  172.       THK(M1)=XY
  173.       GO TO 42
  174. C
  175. C     LOOK FOR THE SECTIN SHEAD WORD 'MATERIAL'
  176. C
  177. 50    REWIND 19
  178. 51    READ(19,900) TEST
  179.       IF(TEST.NE.TEST3) GO TO 51
  180. C
  181. C     BEGIN TO READ IN MATERIAL INFORMATION
  182. C
  183. 52    READ(19,930) M1,AA
  184.       IF(M1.LE.0) GO TO 60
  185.       IF(AA(5).GT.0.0) RHO(M1)=AA(5)
  186.       IF(AA(6).GT.0.0) CP(M1)=AA(6)
  187.       IF(AA(7).GT.0.0) TK(M1)=AA(7)
  188.       CON(18)=RHO(M1)
  189.       CON(19)=CP(M1)
  190.       CON(20)=TK(M1)
  191.       GO TO 52
  192. 60    REWIND 19
  193.       DO 61 JN=1,8
  194. 61    IX(JN)=0
  195.       M1=0
  196.       M2=0
  197. C
  198. C     SEARCH THE SECTION HEAD WORD 'ELEMENTS'
  199. C
  200. 62    READ(19,900) TEST
  201.       IF(TEST.NE.TEST1) GO TO 62
  202. 65    DO 63 JJ=1,8
  203. 63    LX(JJ)=IX(JJ)
  204.       LM1=M1
  205. C
  206. C     BEGIN TO READ ELEMENT INFORMATION AND CALCULATE NETWORD PROPERTIES
  207. C
  208. 64    READ(19,910) M1,M2,(IX(J),J=1,8),L1,L2,L3,L4,K1,K2
  209.       WRITE(2,910) M1,M2,(IX(J),J=1,8),L1,L2,L3,L4,K1,K2
  210.       IF(K1*K2.NE.0.AND.K2.LT.K1) GO TO 998
  211.       IF(M1.LE.0) GO TO 999
  212.       IF(LM1.EQ.0) LM1=M1
  213.       IF(M1-LM1.LE.0) GO TO 65
  214.       GO TO (100,100,200,200,300,200,998,200,100,300,300,300),M2
  215. C
  216. C     IF M2 EQUALS 1,2,9, THEN MAKE ONE DIMENSINAL CALCULATION FOR THE
  217. C     TRUSS.BEAM AND CURVED ELBOW ELEMENT.
  218. C
  219. 100   DO 120 II=LM1,M1  
  220. C     IF(II.EQ.1) GO TO 199
  221. C     IF(LX(1).EQ.0.AND.LX(2).EQ.0) GO TO 199
  222.       IX(1)=LX(1)+(II-LM1)*K1
  223.       IX(2)=LX(2)+(II-LM1)*K1
  224. C     WRITE(2,198) IX(1),IX(2)
  225. C198   FORMAT(1X,2I5)
  226.       J=IX(1)
  227.       KJ=IX(2)
  228.       DX1=A(J)-A(KJ)
  229.       DX2=B(J)-B(KJ)
  230.       DX3=D(J)-D(KJ)
  231.       DL=SQRT(DX1*DX1+DX2*DX2+DX3*DX3)
  232.       IF(M2.EQ.1) GO TO 105
  233.       ASEC=BSECA(L1)
  234.       MNUM=MBEAM(L1)
  235.       GO TO 106
  236. 105   ASEC=TSECA(L1)
  237.       MNUM=MTRUS(L1)
  238. 106   CALL SEARCH(MAXNOY,INP,J,KJ,M)
  239.       Y(INP)=Y(INP)+TK(MNUM)*ASEC/DL
  240.       NY1(INP)=M*J+KJ
  241.       XY=0.5*DL*ASEC*RHO(MNUM)*CP(MNUM)
  242.       C(J)=C(J)+XY
  243.       C(KJ)=C(KJ)+XY
  244. C110   LX(1)=IX(1)
  245. C      LX(2)=IX(2)
  246. 120   CONTINUE
  247.       GO TO 65
  248. C
  249. C     IF M2 EQUALS 3,4,6,8. THEN MAKE TWO DIMENSIONAL CALCULATIONS FOR
  250. C     PLANE STRESS, AXISYMMETRIC, PLATE/SHELL AND PLANE STRAIN ELEMENTS
  251. C
  252. 200   IF(K2.EQ.0) IMAX=1
  253.       IF(K2.EQ.0) JMAX=M1-LM1+1
  254.       IF(K2.NE.0) JMAX=K2/K1-1
  255.       IF(K2.NE.0) IMAX=(M1-LM1+1)/JMAX
  256.       IF(K1.EQ.0) K1=IX(1)-LX(1)
  257.       DO 220 I=1,IMAX
  258.       DO 220 J=1,JMAX
  259. C     IF(I.EQ.1.AND.J.EQ.1) GO TO 220
  260.       K=K2*(I-1)+K1*(J-1)
  261.       IN=LX(1)+K
  262.       JN=LX(2)+K
  263.       KN=LX(3)+K
  264.       LN=LX(4)+K
  265. C      WRITE(2,133) IN,JN,KN,LN
  266. C133      FORMAT(1X,4I5)
  267. C     CALL TRIGON(AREA,ANGL,IN,JN,KN)
  268. C     XY=THK(L2)*RHO(L1)*CP(L1)
  269. C     KN=LX(3)+K
  270. C     LN=LX(4)+K
  271.       CALL TRIGON(AREA,ANGL,IN,JN,KN)
  272.       XY=THK(L2)*RHO(L1)*CP(L1)
  273.       C(IN)=C(IN)+AREA(1)*XY
  274.       C(JN)=C(JN)+AREA(2)*XY
  275.       C(KN)=C(KN)+AREA(3)*XY
  276.       CALL SEARCH(MAXNOY,INP,IN,JN,M)
  277.       XY= 0.5*TK(L1)*THK(L2)
  278.       X=ANGL(3)
  279.       Y(INP)=Y(INP)+XY*COT(X)
  280.       NY1(INP)=M*IN+JN
  281.       CALL SEARCH(MAXNOY,INP,JN,KN,M)
  282.       X=ANGL(1)
  283.       Y(INP)=Y(INP)+XY*COT(X)
  284.       NY1(INP)=M*JN+KN
  285.       CALL SEARCH(MAXNOY,INP,KN,IN,M)
  286.       X=ANGL(2)
  287.       Y(INP)=Y(INP)+XY*COT(X)
  288.       NY1(INP)=M*KN+IN
  289.       IF(LN.EQ.0) GO TO 220
  290.       CALL TRIGON(AREA,ANGL,IN,KN,LN)
  291.       XY=THK(L2)*RHO(L1)*CP(L1)
  292.       C(IN)=C(IN)+AREA(1)*XY
  293.       C(KN)=C(KN)+AREA(2)*XY
  294.       C(LN)=C(LN)+AREA(3)*XY
  295.       CALL SEARCH(MAXNOY,INP,IN,KN,M)
  296.       XY=0.5*TK(L1)*THK(L2)
  297.       X=ANGL(3)
  298.       Y(INP)=Y(INP)+XY*COT(X)
  299.          NY1(INP)=M*IN+KN
  300.       CALL SEARCH(MAXNOY,INP,KN,LN,M)
  301.       X=ANGL(1)
  302.       Y(INP)=Y(INP)+XY*COT(X)
  303.       NY1(INP)=M*KN+LN
  304.       CALL SEARCH(MAXNOY,INP,LN,IN,M)
  305.       X=ANGL(2)
  306.       Y(INP)=Y(INP)+XY*COT(X)
  307.       NY1(INP)=M*LN+IN
  308. 220   CONTINUE
  309.       GO TO 65
  310. C
  311. C     IF M2 EQUALS 5,10, THEN MAKE THREE DIMENSIONAL CALCULATIONS FOR
  312. C     SOLID ELEMENTS
  313. C
  314. 300   IF(K2.EQ.0) IMAX=1
  315.       IF(K2.EQ.0) JMAX=M1-LM1+1
  316.       IF(K2.NE.0) JMAX=K2/K1-1
  317.       IF(K2.NE.0) IMAX=(M1-LM1+1)/JMAX
  318.       IF(K1.EQ.0) K1=IX(1)-LX(1)
  319.       DO 320 I=1,IMAX
  320.       DO 320 J=1,JMAX
  321. C     IF(I.EQ.1.AND.J.EQ.1) GO TO 320
  322.       K=K2*(I-1)+K1*(J-1)
  323.       IN1=LX(1)+K
  324.       JN1=LX(2)+K
  325.       KN1=LX(3)+K
  326.       LN1=LX(4)+K
  327.       IN2=LX(5)+K
  328.       JN2=LX(6)+K
  329.       KN2=LX(7)+K
  330.       LN2=LX(8)+K
  331.       DL1=SQRT((A(IN1)-A(IN2))**2+(B(IN1)-B(IN2))**2+(D(IN1)-D(IN2))**2)
  332.       DL2=SQRT((A(JN1)-A(JN2))**2+(B(JN1)-B(JN2))**2+(D(JN1)-D(JN2))**2)
  333.       DL3=SQRT((A(KN1)-A(KN2))**2+(B(KN1)-B(KN2))**2+(D(KN1)-D(KN2))**2)
  334.       DL4=SQRT((A(LN1)-A(LN2))**2+(B(LN1)-B(LN2))**2+(D(LN1)-D(LN2))**2)
  335.       XY=RHO(L1)*CP(L1)
  336.       XYZ=0.5*TK(L1)
  337.       CALL TRIGON(AREA,ANGL,IN1,JN1,KN1)
  338.       C(IN1)=C(IN1)+AREA(1)*XY*0.5*DL1
  339.       C(JN1)=C(JN1)+AREA(2)*XY*0.5*DL2
  340.       C(KN1)=C(KN1)+AREA(3)*XY*0.5*DL3
  341.       CALL SEARCH(MAXNOY,INP,IN1,JN1,M)
  342.       X=ANGL(3)
  343.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
  344.       NY1(INP)=M*IN1+JN1
  345.       CALL SEARCH(MAXNOY,INP,JN1,KN1,M)
  346.       X=ANGL(1)
  347.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL2
  348.       NY1(INP)=M*JN1+KN1
  349.       CALL SEARCH(MAXNOY,INP,KN1,IN1,M)
  350.       X=ANGL(2)
  351.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
  352.       NY1(INP)=M*KN1+IN1
  353.       CALL SEARCH(MAXNOY,INP,IN1,IN2,M)
  354.       Y(INP)=Y(INP)+TK(L1)*AREA(1)/DL1
  355.       NY1(INP)=M*IN1+IN2
  356.       CALL SEARCH(MAXNOY,INP,JN1,JN2,M)
  357.       Y(INP)=Y(INP)+TK(L1)*AREA(2)/DL2
  358.       NY1(INP)=M*JN1+JN2
  359.       CALL SEARCH(MAXNOY,INP,KN1,KN2,M)
  360.       Y(INP)=Y(INP)+TK(L1)*AREA(3)/DL3
  361.       NY1(INP)=M*KN1+KN2
  362.       CALL TRIGON(AREA,ANGL,IN1,KN1,LN1)
  363.       C(IN1)=C(IN1)+AREA(1)*XY*0.5*DL1
  364.       C(KN1)=C(KN1)+AREA(2)*XY*0.5*DL3
  365.       C(LN1)=C(LN1)+AREA(3)*XY*0.5*DL4
  366.       CALL SEARCH(MAXNOY,INP,IN1,KN1,M)
  367.       X=ANGL(3)
  368.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
  369.       NY1(INP)=M*IN1+KN1
  370.       CALL SEARCH(MAXNOY,INP,KN1,LN1,M)
  371.       X=ANGL(1)
  372.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
  373.       NY1(INP)=M*KN1+LN1
  374.       CALL SEARCH(MAXNOY,INP,LN1,IN1,M)
  375.       X=ANGL(2)
  376.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL4
  377.       NY1(INP)=M*LN1+IN1
  378.       CALL SEARCH(MAXNOY,INP,IN1,IN2,M)
  379.       Y(INP)=Y(INP)+TK(L1)*AREA(1)/DL1
  380.       NY1(INP)=M*IN1+IN2
  381.       CALL SEARCH(MAXNOY,INP,KN1,KN2,M)
  382.       Y(INP)=Y(INP)+TK(L1)*AREA(2)/DL3
  383.       NY1(INP)=M*KN1+KN2
  384.       CALL SEARCH(MAXNOY,INP,LN1,LN2,M)
  385.       Y(INP)=Y(INP)+TK(L1)*AREA(3)/DL4
  386.       NY1(INP)=M*LN1+LN2
  387.       CALL TRIGON(AREA,ANGL,IN2,JN2,KN2)
  388.       C(IN2)=C(IN2)+AREA(1)*XY*0.5*DL1
  389.       C(JN2)=C(JN2)+AREA(2)*XY*0.5*DL2
  390.       C(KN2)=C(KN2)+AREA(3)*XY*0.5*DL3
  391.       CALL SEARCH(MAXNOY,INP,IN2,JN2,M)
  392.       X=ANGL(3)
  393.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
  394.       NY1(INP)=M*IN2+JN2
  395.       CALL SEARCH(MAXNOY,INP,JN2,KN2,M)
  396.       X=ANGL(1)
  397.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL2
  398.       NY1(INP)=M*JN2+KN2
  399.       CALL SEARCH(MAXNOY,INP,KN2,IN2,M)
  400.       X=ANGL(2)
  401.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
  402.       NY1(INP)=M*KN2+IN2
  403.       CALL TRIGON(AREA,ANGL,IN2,KN2,LN2)
  404.       C(IN2)=C(IN2)+AREA(1)*XY*0.5*DL3
  405.       C(KN2)=C(KN2)+AREA(2)*XY*0.5*DL3
  406.       C(LN2)=C(LN2)+AREA(3)*XY*0.5*DL4
  407.       CALL SEARCH(MAXNOY,INP,IN2,KN2,M)
  408.       X=ANGL(3)
  409.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
  410.       NY1(INP)=M*IN2+KN2
  411.       CALL SEARCH(MAXNOY,INP,LN2,IN2,M)
  412.       X=ANGL(1)
  413.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
  414.       NY1(INP)=M*LN2+IN2
  415.       CALL SEARCH(MAXNOY,INP,KN2,LN2,M)
  416.       X=ANGL(2)
  417.       Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL4
  418.       NY1(INP)=M*KN2+LN2
  419. 320   CONTINUE
  420.       GO TO 65
  421. 900   FORMAT(A8)
  422. 920   FORMAT(I10,3F10.0,I5,4X,A1,I5,F10.0)
  423. 970   FORMAT(2I10,F10.0)
  424. 930   FORMAT(I10,7F10.0)
  425. 910   FORMAT(16I5)
  426. 960   FORMAT(1X,'ERROR IN CONY',2X, 'K2 MUST BE GREATER THAN K1')
  427. 998   WRITE(6,960)
  428.       CALL EXIT
  429. 999   RETURN
  430.       END
  431. C     SUBROUTINE TRIGON IS USED TO CALCULATE THE AERA AND ANGLES OF A
  432. C     TRIANGLE
  433. C
  434.       SUBROUTINE TRIGON(AREA,ANGL,I,J,K)
  435.       DIMENSION  ZZ(16000),A(1000),B(1000),D(1000),AREA(3),ANGL(3)
  436.       COMMON ZZ
  437.       EQUIVALENCE (ZZ(7002),A(1)),(ZZ(8002),B(1)),(ZZ(6002),D(1))
  438.       AX=A(J)-A(K)
  439.       AY=B(J)-B(K)
  440.       AZ=D(J)-D(K)
  441.       AL=SQRT(AX*AX+AY*AY+AZ*AZ)
  442.       BX=A(K)-A(I)
  443.       BY=B(K)-B(I)
  444.       BZ=D(K)-D(I)
  445.       BL=SQRT(BX*BX+BY*BY+BZ*BZ)
  446.       CX=A(I)-A(J)
  447.       CY=B(I)-B(J)
  448.       CZ=D(I)-D(J)
  449.       CL=SQRT(CX*CX+CY*CY+CZ*CZ)
  450.       X=((BX*CX+BY*CY+BZ*CZ)/(BL*CL))
  451.       IF(X.LT.0.0) X=-X
  452.       ANGL(1)=ACOS(X)
  453.       X=(CX*AX+CY*AY+CZ*AZ)/(CL*AL)
  454.       IF(X.LT.0.0) X=-X
  455.       ANGL(2)=ACOS(X)
  456.       X=(AX*BX+AY*BY+AZ*BZ)/(AL*BL)
  457.       IF(X.LT.0.0) X=-X
  458.       ANGL(3)=ACOS(X)
  459.       X=ANGL(2)
  460.       AREA(1)=0.125*BL*BL*COT(X)
  461.       X=ANGL(3)
  462.       AREA(2)=0.125*CL*CL*COT(X)
  463.       X=ANGL(1)
  464.       AREA(3)=0.125*AL*AL*COT(X)
  465.       X=AREA(1)
  466.       AREA(1)=AREA(1)+AREA(2)
  467.       AREA(2)=AREA(2)+AREA(3)
  468.       AREA(3)=AREA(3)+X
  469.       RETURN
  470.       END
  471.