C     SUBROUTINE CONY  CONECT SAP6 WITH TAP6  IT IS USED TO GENERATION  
C     TEMPRATURE TO SAP6 AND THROUGH SAP6 CALCULAT THERMAL STREE AND
C     DISPLACEMENT. CONY CAN READ SAP6 DATA CARD AND TRANFORM TO TAP6 DIF
C     ERECE NETWORK. WE GET THIS VERSION SUB CONY HAS MANY ERROR SO CAN'T
C     BE USED.
C      THIS SUB IS MODIFIED BY T.U.P  CAD OFFICE IN BEIJING, 1986
C
      SUBROUTINE CONY
      INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000),
     X        NXP(1000),KEY(1000),IZ(1000),IPR(500),IPH(100),CON(20),
     X        NTHEAD(100),NTTAIL(100),NFLAG(1000),MBEAM(100),MTRUS(100)
      DIMENSION ZZ(16000),DYTCQ(11030),Y(3000),T(1000),C(1000),Q(1000),
     X          GAIN(500),A(1000),B(1000),D(1000),RHO(50),CP(50),TK(50),
     X          TPROG(2),AA(7),IX(8),LX(8),XX(3),XO(3),KS(4),THK(100),
     X          BSECA(100),TSECA(100),AREA(3),ANGL(3)
      REAL*8 TEST,TEST1,TEST2,TEST3,TEST4,TEST5,TEST6,TEST7,TEST8,TEST9
      COMMON /I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IRP,IPH,NTHEAD,
     X        NTTAIL,NFLAG,IZ
      COMMON /LAB1/PLOT(100,10),IPLOT(10)
      COMMON /LAB2/DATE(2)
      COMMON /LAB3/ICASE,XLIM(2),YLIM(2),CT(5),XLAB(5),YLAB(5)
      COMMON /LAB4/SPROG(2)
      COMMON /LAB5/TITLE
      COMMON ZZ
      EQUIVALENCE (ZZ(1),TIME),(ZZ(2),Y(1),DYTCQ(1)),
     X    (ZZ(3002),T(1)),(ZZ(4002),C(1)),(ZZ(11061),MAX1T)
     X           ,(ZZ(5002),Q(1)),(ZZ(6002),D(1)),(ZZ(11060),MIN1T),
     X    (ZZ(7002),A(1)),(ZZ(8002),B(1)),(ZZ(11077),IMIN),
     X    (ZZ(11078),IMAX),(ZZ(11085),MAXNOT),(ZZ(11086),MAXNOY),
     X    (ZZ(11002),CON(1))
      DATA TEST1/8HELEMENTS/
      DATA TEST2/8HNODE-INP/
      DATA TEST3/8HMATERIAL/
      DATA TEST4/8HEXECUTE-/
      DATA TEST5/8HREADDATA/
      DATA TEST6/8HTHIK-TYP/
      DATA TEST7/8HTRUS-TYP/
      DATA TEST8/8HBEAM-TYP/
      DATA TEST9/8HTEMPDIST/
      DATA KS/1H ,1HR,1HC,1HS/
      REWIND 19
      M=1000
      PI2=(3.14159/180.0)**2
C
C     SEAARCH THE SECTION HEAD WORD 'NODE-INP'
C
10    READ(19,900) TEST
      IF(TEST.NE.TEST2) GO TO 10
      DO 2  J=1,8
      IX(J)=0
2     CONTINUE
          DO 134 JJ=1,3
134       XX(JJ)=0.0
11    DO 5 J=1,3
      XO(J)=XX(J)
5     LX(J)=IX(J)
C
C     BEGIN TO READ NODE INFORMATION
C
      READ(19,920) IX(1),XX(1),XX(2),XX(3),IX(2),KT,IX(3),XY
      IF(IX(1).LE.0) GO TO 20
      L2=IX(1)
      MAX1T=MAX0(MAX1T,L2)
      MIN1T=MIN0(MIN1T,L2)
      MAXNOT=MAX1T
      IF(KT.EQ.KS(3)) GO TO 6
      IF(KT.EQ.KS(4)) GO TO 7
      NFLAG(L2)=2
      A(L2)=XX(1)
      B(L2)=XX(2)
      D(L2)=XX(3)
      GO TO 8
6     NFLAG(L2)=3
      PHI=XX(2)
      A(L2)=XX(1)*COS(PHI)
      B(L2)=XX(1)*SIN(PHI)
      D(L2)=XX(3)
      GO TO 8
7     NFLAG(L2)=4
      PHI=XX(2)
      THETA=XX(3)
      A(L2)=XX(1)*SIN(THETA)*COS(PHI)
      B(L2)=XX(1)*SIN(THETA)*SIN(PHI)
      D(L2)=XX(1)*COS(THETA)
8     IF(IX(3).LE.0) GO TO 11
      L1=LX(1)
      L3=IX(3)
      DN=0.0
      X=0.01*XY
      IF(XY.LE.0.0) X=1.0
      LDX=L2-L1
      K=0
      DO 12 J=1,LDX,L3
      K=K+1
12    DN=DN+X**(K-1)
      DX1=(XX(1)-XO(1))/DN
      DX2=(XX(2)-XO(2))/DN
      DX3=(XX(3)-XO(3))/DN
      K=0
      DO 15 J=L1,L2,L3
      K=K+1
      IF(J.EQ.1) GO TO 15
      XX(1)=XO(1)+DX1*(K-1)*X**(J-L1)
      XX(2)=XO(2)+DX2*(K-1)*X**(J-L1)
      XX(3)=XO(3)+DX3*(K-1)*X**(J-L1)
132   IF(KT.EQ.KS(3)) GO TO 16
      IF(KT.EQ.KS(4)) GO TO 17
      NFLAG(J)=2
      A(J)=XX(1)
      B(J)=XX(2)
      D(J)=XX(3)
      GO TO 15
16    NFLAG(J)=3
      PHI=XX(2)
      A(J)=XX(1)*COS(PHI)
      B(J)=XX(1)*SIN(PHI)
      D(J)=XX(3)
      GO TO 15
17    NFLAG(J)=4
      PHI=XX(2)
      THETA=XX(3)
      A(J)=XX(1)*SIN(THETA)*COS(PHI)
      B(J)=XX(1)*SIN(THETA)*SIN(PHI)
      D(J)=XX(1)*COS(THETA)
C8    XO(1)=XX(1)
C     XO(2)=XX(2)
C     XO(3)=XX(3)
15    CONTINUE
      GO TO 11
20    REWIND 19
C
C     SEARCH THE SECTION HEAD WORD 'TRUS-TYP'
C
21    READ(19,900) TEST
      IF(TEST.EQ.TEST4.OR.TEST.EQ.TEST5) GO TO 30
      IF(TEST.NE.TEST7) GO TO 21
C
C     BEGIN TO READ TRUS-TYP PROERTIES
C
22    READ(19,970) M1,M2,XY
      IF(M1.LE.0) GO TO 30
      TSECA(M1)=XY
      MTRUS(M1)=M2
      GO TO 22
30    REWIND 19
C
C     SEARCH THE SECTION HEAD WORD 'BEAM-TYP'
C
31    READ(19,900) TEST
      IF(TEST.EQ.TEST4.OR.TEST.EQ.TEST5) GO TO 40
      IF(TEST.NE.TEST8) GO TO 31
C
C     BEGIN TO READ BEAM-TYP PROPERTIES
C
32    READ(19,970) M1,M2,XY
      IF(M1.LE.0) GO TO 40
      BSECA(M1)=XY
      MBEAM(M1)=M2
      GO TO 32
40    REWIND 19
C
C     SEARCH THE SECTION HEAD WORD 'THIK-TYP'
C
41    READ(19,900) TEST
      IF(TEST.EQ.TEST4.OR.TEST.EQ.TEST5) GO TO 50
      IF(TEST.NE.TEST6) GO TO 41
C
C     BEGIN TO READ THE THICKNESS OF PLATE-SHELL EMEMENTS
C
42    READ(19,930) M1,XY
      IF(M1.LE.0) GO TO 50
      THK(M1)=XY
      GO TO 42
C
C     LOOK FOR THE SECTIN SHEAD WORD 'MATERIAL'
C
50    REWIND 19
51    READ(19,900) TEST
      IF(TEST.NE.TEST3) GO TO 51
C
C     BEGIN TO READ IN MATERIAL INFORMATION
C
52    READ(19,930) M1,AA
      IF(M1.LE.0) GO TO 60
      IF(AA(5).GT.0.0) RHO(M1)=AA(5)
      IF(AA(6).GT.0.0) CP(M1)=AA(6)
      IF(AA(7).GT.0.0) TK(M1)=AA(7)
      CON(18)=RHO(M1)
      CON(19)=CP(M1)
      CON(20)=TK(M1)
      GO TO 52
60    REWIND 19
      DO 61 JN=1,8
61    IX(JN)=0
      M1=0
      M2=0
C
C     SEARCH THE SECTION HEAD WORD 'ELEMENTS'
C
62    READ(19,900) TEST
      IF(TEST.NE.TEST1) GO TO 62
65    DO 63 JJ=1,8
63    LX(JJ)=IX(JJ)
      LM1=M1
C
C     BEGIN TO READ ELEMENT INFORMATION AND CALCULATE NETWORD PROPERTIES
C
64    READ(19,910) M1,M2,(IX(J),J=1,8),L1,L2,L3,L4,K1,K2
      WRITE(2,910) M1,M2,(IX(J),J=1,8),L1,L2,L3,L4,K1,K2
      IF(K1*K2.NE.0.AND.K2.LT.K1) GO TO 998
      IF(M1.LE.0) GO TO 999
      IF(LM1.EQ.0) LM1=M1
      IF(M1-LM1.LE.0) GO TO 65
      GO TO (100,100,200,200,300,200,998,200,100,300,300,300),M2
C
C     IF M2 EQUALS 1,2,9, THEN MAKE ONE DIMENSINAL CALCULATION FOR THE
C     TRUSS.BEAM AND CURVED ELBOW ELEMENT.
C
100   DO 120 II=LM1,M1  
C     IF(II.EQ.1) GO TO 199
C     IF(LX(1).EQ.0.AND.LX(2).EQ.0) GO TO 199
      IX(1)=LX(1)+(II-LM1)*K1
      IX(2)=LX(2)+(II-LM1)*K1
C     WRITE(2,198) IX(1),IX(2)
C198   FORMAT(1X,2I5)
      J=IX(1)
      KJ=IX(2)
      DX1=A(J)-A(KJ)
      DX2=B(J)-B(KJ)
      DX3=D(J)-D(KJ)
      DL=SQRT(DX1*DX1+DX2*DX2+DX3*DX3)
      IF(M2.EQ.1) GO TO 105
      ASEC=BSECA(L1)
      MNUM=MBEAM(L1)
      GO TO 106
105   ASEC=TSECA(L1)
      MNUM=MTRUS(L1)
106   CALL SEARCH(MAXNOY,INP,J,KJ,M)
      Y(INP)=Y(INP)+TK(MNUM)*ASEC/DL
      NY1(INP)=M*J+KJ
      XY=0.5*DL*ASEC*RHO(MNUM)*CP(MNUM)
      C(J)=C(J)+XY
      C(KJ)=C(KJ)+XY
C110   LX(1)=IX(1)
C      LX(2)=IX(2)
120   CONTINUE
      GO TO 65
C
C     IF M2 EQUALS 3,4,6,8. THEN MAKE TWO DIMENSIONAL CALCULATIONS FOR
C     PLANE STRESS, AXISYMMETRIC, PLATE/SHELL AND PLANE STRAIN ELEMENTS
C
200   IF(K2.EQ.0) IMAX=1
      IF(K2.EQ.0) JMAX=M1-LM1+1
      IF(K2.NE.0) JMAX=K2/K1-1
      IF(K2.NE.0) IMAX=(M1-LM1+1)/JMAX
      IF(K1.EQ.0) K1=IX(1)-LX(1)
      DO 220 I=1,IMAX
      DO 220 J=1,JMAX
C     IF(I.EQ.1.AND.J.EQ.1) GO TO 220
      K=K2*(I-1)+K1*(J-1)
      IN=LX(1)+K
      JN=LX(2)+K
      KN=LX(3)+K
      LN=LX(4)+K
C      WRITE(2,133) IN,JN,KN,LN
C133      FORMAT(1X,4I5)
C     CALL TRIGON(AREA,ANGL,IN,JN,KN)
C     XY=THK(L2)*RHO(L1)*CP(L1)
C     KN=LX(3)+K
C     LN=LX(4)+K
      CALL TRIGON(AREA,ANGL,IN,JN,KN)
      XY=THK(L2)*RHO(L1)*CP(L1)
      C(IN)=C(IN)+AREA(1)*XY
      C(JN)=C(JN)+AREA(2)*XY
      C(KN)=C(KN)+AREA(3)*XY
      CALL SEARCH(MAXNOY,INP,IN,JN,M)
      XY= 0.5*TK(L1)*THK(L2)
      X=ANGL(3)
      Y(INP)=Y(INP)+XY*COT(X)
      NY1(INP)=M*IN+JN
      CALL SEARCH(MAXNOY,INP,JN,KN,M)
      X=ANGL(1)
      Y(INP)=Y(INP)+XY*COT(X)
      NY1(INP)=M*JN+KN
      CALL SEARCH(MAXNOY,INP,KN,IN,M)
      X=ANGL(2)
      Y(INP)=Y(INP)+XY*COT(X)
      NY1(INP)=M*KN+IN
      IF(LN.EQ.0) GO TO 220
      CALL TRIGON(AREA,ANGL,IN,KN,LN)
      XY=THK(L2)*RHO(L1)*CP(L1)
      C(IN)=C(IN)+AREA(1)*XY
      C(KN)=C(KN)+AREA(2)*XY
      C(LN)=C(LN)+AREA(3)*XY
      CALL SEARCH(MAXNOY,INP,IN,KN,M)
      XY=0.5*TK(L1)*THK(L2)
      X=ANGL(3)
      Y(INP)=Y(INP)+XY*COT(X)
         NY1(INP)=M*IN+KN
      CALL SEARCH(MAXNOY,INP,KN,LN,M)
      X=ANGL(1)
      Y(INP)=Y(INP)+XY*COT(X)
      NY1(INP)=M*KN+LN
      CALL SEARCH(MAXNOY,INP,LN,IN,M)
      X=ANGL(2)
      Y(INP)=Y(INP)+XY*COT(X)
      NY1(INP)=M*LN+IN
220   CONTINUE
      GO TO 65
C
C     IF M2 EQUALS 5,10, THEN MAKE THREE DIMENSIONAL CALCULATIONS FOR
C     SOLID ELEMENTS
C
300   IF(K2.EQ.0) IMAX=1
      IF(K2.EQ.0) JMAX=M1-LM1+1
      IF(K2.NE.0) JMAX=K2/K1-1
      IF(K2.NE.0) IMAX=(M1-LM1+1)/JMAX
      IF(K1.EQ.0) K1=IX(1)-LX(1)
      DO 320 I=1,IMAX
      DO 320 J=1,JMAX
C     IF(I.EQ.1.AND.J.EQ.1) GO TO 320
      K=K2*(I-1)+K1*(J-1)
      IN1=LX(1)+K
      JN1=LX(2)+K
      KN1=LX(3)+K
      LN1=LX(4)+K
      IN2=LX(5)+K
      JN2=LX(6)+K
      KN2=LX(7)+K
      LN2=LX(8)+K
      DL1=SQRT((A(IN1)-A(IN2))**2+(B(IN1)-B(IN2))**2+(D(IN1)-D(IN2))**2)
      DL2=SQRT((A(JN1)-A(JN2))**2+(B(JN1)-B(JN2))**2+(D(JN1)-D(JN2))**2)
      DL3=SQRT((A(KN1)-A(KN2))**2+(B(KN1)-B(KN2))**2+(D(KN1)-D(KN2))**2)
      DL4=SQRT((A(LN1)-A(LN2))**2+(B(LN1)-B(LN2))**2+(D(LN1)-D(LN2))**2)
      XY=RHO(L1)*CP(L1)
      XYZ=0.5*TK(L1)
      CALL TRIGON(AREA,ANGL,IN1,JN1,KN1)
      C(IN1)=C(IN1)+AREA(1)*XY*0.5*DL1
      C(JN1)=C(JN1)+AREA(2)*XY*0.5*DL2
      C(KN1)=C(KN1)+AREA(3)*XY*0.5*DL3
      CALL SEARCH(MAXNOY,INP,IN1,JN1,M)
      X=ANGL(3)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
      NY1(INP)=M*IN1+JN1
      CALL SEARCH(MAXNOY,INP,JN1,KN1,M)
      X=ANGL(1)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL2
      NY1(INP)=M*JN1+KN1
      CALL SEARCH(MAXNOY,INP,KN1,IN1,M)
      X=ANGL(2)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
      NY1(INP)=M*KN1+IN1
      CALL SEARCH(MAXNOY,INP,IN1,IN2,M)
      Y(INP)=Y(INP)+TK(L1)*AREA(1)/DL1
      NY1(INP)=M*IN1+IN2
      CALL SEARCH(MAXNOY,INP,JN1,JN2,M)
      Y(INP)=Y(INP)+TK(L1)*AREA(2)/DL2
      NY1(INP)=M*JN1+JN2
      CALL SEARCH(MAXNOY,INP,KN1,KN2,M)
      Y(INP)=Y(INP)+TK(L1)*AREA(3)/DL3
      NY1(INP)=M*KN1+KN2
      CALL TRIGON(AREA,ANGL,IN1,KN1,LN1)
      C(IN1)=C(IN1)+AREA(1)*XY*0.5*DL1
      C(KN1)=C(KN1)+AREA(2)*XY*0.5*DL3
      C(LN1)=C(LN1)+AREA(3)*XY*0.5*DL4
      CALL SEARCH(MAXNOY,INP,IN1,KN1,M)
      X=ANGL(3)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
      NY1(INP)=M*IN1+KN1
      CALL SEARCH(MAXNOY,INP,KN1,LN1,M)
      X=ANGL(1)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
      NY1(INP)=M*KN1+LN1
      CALL SEARCH(MAXNOY,INP,LN1,IN1,M)
      X=ANGL(2)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL4
      NY1(INP)=M*LN1+IN1
      CALL SEARCH(MAXNOY,INP,IN1,IN2,M)
      Y(INP)=Y(INP)+TK(L1)*AREA(1)/DL1
      NY1(INP)=M*IN1+IN2
      CALL SEARCH(MAXNOY,INP,KN1,KN2,M)
      Y(INP)=Y(INP)+TK(L1)*AREA(2)/DL3
      NY1(INP)=M*KN1+KN2
      CALL SEARCH(MAXNOY,INP,LN1,LN2,M)
      Y(INP)=Y(INP)+TK(L1)*AREA(3)/DL4
      NY1(INP)=M*LN1+LN2
      CALL TRIGON(AREA,ANGL,IN2,JN2,KN2)
      C(IN2)=C(IN2)+AREA(1)*XY*0.5*DL1
      C(JN2)=C(JN2)+AREA(2)*XY*0.5*DL2
      C(KN2)=C(KN2)+AREA(3)*XY*0.5*DL3
      CALL SEARCH(MAXNOY,INP,IN2,JN2,M)
      X=ANGL(3)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
      NY1(INP)=M*IN2+JN2
      CALL SEARCH(MAXNOY,INP,JN2,KN2,M)
      X=ANGL(1)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL2
      NY1(INP)=M*JN2+KN2
      CALL SEARCH(MAXNOY,INP,KN2,IN2,M)
      X=ANGL(2)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
      NY1(INP)=M*KN2+IN2
      CALL TRIGON(AREA,ANGL,IN2,KN2,LN2)
      C(IN2)=C(IN2)+AREA(1)*XY*0.5*DL3
      C(KN2)=C(KN2)+AREA(2)*XY*0.5*DL3
      C(LN2)=C(LN2)+AREA(3)*XY*0.5*DL4
      CALL SEARCH(MAXNOY,INP,IN2,KN2,M)
      X=ANGL(3)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL1
      NY1(INP)=M*IN2+KN2
      CALL SEARCH(MAXNOY,INP,LN2,IN2,M)
      X=ANGL(1)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL3
      NY1(INP)=M*LN2+IN2
      CALL SEARCH(MAXNOY,INP,KN2,LN2,M)
      X=ANGL(2)
      Y(INP)=Y(INP)+XYZ*COT(X)*0.5*DL4
      NY1(INP)=M*KN2+LN2
320   CONTINUE
      GO TO 65
900   FORMAT(A8)
920   FORMAT(I10,3F10.0,I5,4X,A1,I5,F10.0)
970   FORMAT(2I10,F10.0)
930   FORMAT(I10,7F10.0)
910   FORMAT(16I5)
960   FORMAT(1X,'ERROR IN CONY',2X, 'K2 MUST BE GREATER THAN K1')
998   WRITE(6,960)
      CALL EXIT
999   RETURN
      END
C     SUBROUTINE TRIGON IS USED TO CALCULATE THE AERA AND ANGLES OF A
C     TRIANGLE
C
      SUBROUTINE TRIGON(AREA,ANGL,I,J,K)
      DIMENSION  ZZ(16000),A(1000),B(1000),D(1000),AREA(3),ANGL(3)
      COMMON ZZ
      EQUIVALENCE (ZZ(7002),A(1)),(ZZ(8002),B(1)),(ZZ(6002),D(1))
      AX=A(J)-A(K)
      AY=B(J)-B(K)
      AZ=D(J)-D(K)
      AL=SQRT(AX*AX+AY*AY+AZ*AZ)
      BX=A(K)-A(I)
      BY=B(K)-B(I)
      BZ=D(K)-D(I)
      BL=SQRT(BX*BX+BY*BY+BZ*BZ)
      CX=A(I)-A(J)
      CY=B(I)-B(J)
      CZ=D(I)-D(J)
      CL=SQRT(CX*CX+CY*CY+CZ*CZ)
      X=((BX*CX+BY*CY+BZ*CZ)/(BL*CL))
      IF(X.LT.0.0) X=-X
      ANGL(1)=ACOS(X)
      X=(CX*AX+CY*AY+CZ*AZ)/(CL*AL)
      IF(X.LT.0.0) X=-X
      ANGL(2)=ACOS(X)
      X=(AX*BX+AY*BY+AZ*BZ)/(AL*BL)
      IF(X.LT.0.0) X=-X
      ANGL(3)=ACOS(X)
      X=ANGL(2)
      AREA(1)=0.125*BL*BL*COT(X)
      X=ANGL(3)
      AREA(2)=0.125*CL*CL*COT(X)
      X=ANGL(1)
      AREA(3)=0.125*AL*AL*COT(X)
      X=AREA(1)
      AREA(1)=AREA(1)+AREA(2)
      AREA(2)=AREA(2)+AREA(3)
      AREA(3)=AREA(3)+X
      RETURN
      END
