      PROGRAM WIRE3
C======================================================================
C----- UNEVEN LEMENT<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C---------DZMAXER = MAX ERROR IN Z DURING ITERATIVE PROCEDURE
C                   DZMAXER IS DIVIDED BY DOMAIN LENGTH
C   TX = TENSION IN X-DIRECTION
C        TX IS DIVIDED BY GRAVITY.
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=2, MXE=1000, MXN=MXE+1,NBW=2*(ND-1)+1 )
      PARAMETER ( MXITERA=80, ERRMAX=1.D-4, GRAVITY = 1.D0 )
      DIMENSION NODEX(MXE,ND),X(MXN),Z(MXN),PARA(MXN),
     * STIFF(ND,ND),A(MXN,NBW),RHS(MXN),TS(MXN)
C======================================================================
C (1)-READING OF DATA
      CALL INPUT (NE,DOMAINL,WDENSITY,TXMAX,TXMIN,TXSTEP, Z1,ZN)
      GW = WDENSITY*GRAVITY
C======================================================================
C (2)-FINITE ELEMENT CREATION
      CALL ELEMENT ( DOMAINL, MXE,MXN,ND,NE,NNODE,NODEX,X  )
C======================================================================
C (4)-COMPUTATION OF INITIAL LOCATION OF STRING
C-------INITIAL Z STANDS FOR A SET OF GUESS STRING LOCATION.
      CALL INITIALZ ( MXN,NNODE,GW,TXMAX,X, Z1, ZN, Z )
C======================================================================
      OPEN ( 1,FILE='SOLUTION.FEM', STATUS='UNKNOWN')
      OPEN ( 2,FILE='SUMMARY.FEM', STATUS='UNKNOWN')
      WRITE(1,*) 'X  Z PARABOLIC TENSION TX ITERATIONS'
      WRITE(2,*) 'TX TS-UP TS-DW MAXDEPTH L-OF-WIRE W-OF-WIRE SUM-TZ',
     *           ' %ERROR ITERATIONS'
C
C
      DO TX = TXMAX, TXMIN, -TXSTEP
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      ITERA = 0
      DZMAXER = 2.*ERRMAX
      DO WHILE ( ITERA .LE. MXITERA .AND. DZMAXER .GT. ERRMAX )
C*********************************************************************
      ITERA = ITERA + 1
C==============================================
C (7)-CONSTRUCTION OF FEM-MATRIX EQUATION
      CALL MATRIX ( MXE,MXN,ND,NBW,NE,NNODE, GRAVITY,WDENSITY, 
     * TX,STIFF, NODEX, X, Z, A,RHS )
C==============================================
C (8)-IMPLEMENTATION OF BOUNDARY CONDITION
      CALL FORM ( MXN, NBW, NNODE, A, RHS, Z1, ZN )
C==============================================
C (9)-SOLVE FOR UNKNOWN VARIABLES
      CALL SYSTEM ( MXN, NBW, NNODE, A , RHS )
C==============================================
C(10)-EVALUATION OF NEW STRING LOCATION BASED UPON OLD U AND NEW U.
C-------- UMAX = MAXIMUM ERROR IN U(I) DURING ITERATION
      CALL NEWZ ( MXN,NNODE,DOMAINL,  RHS ,Z, DZMAXER )
      END DO
C-------- END OF ITARETIVE PROCEDURE ---------
C*********************************************************************
      CALL PRT (ND,MXE,MXN,NE,NNODE,NODEX,X,TX,GW,PARA,RHS,TS,
     *                DZMAXER,ITERA)
C------ END OF COMPUTATIONAL PROCEDURE
      END DO
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      CLOSE (2)
      CLOSE (1)
      STOP'NORMAL TERMINATION'
      END
C
C
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE PRT (ND,MXE,MXN,NE,NNODE,NODEX,X,TX,GW,PARA,RHS,TS,
     *                DZMAXER,ITERA)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXN),Z(MXN),PARA(MXN),RHS(MXN),TS(MXN),NODEX(MXE,ND)
      CALL PARAB (MXN,NNODE,X,TX,GW,RHS,PARA)
      CALL CMPDEPTH ( MXN,NNODE,X,RHS,DEPTH,WL )
      CALL TENSION (ND,MXE,MXN,NE,NODEX,X,RHS,TX,TS,SUMTZ)
      WOW = GW*WL
      WRITE(1,*)
      WRITE(1,'(5F13.5,I5)') X(1),RHS(1), PARA(1), TS(1),TX,ITERA
      DO I = 2 , NNODE
      WRITE(1,'(4F13.5)') X(I),RHS(I), PARA(I), TS(I)
      END DO
      PRCNTER = ( WOW - SUMTZ ) / WOW * 100.
      WRITE (2,100) TX,TS(1),TS(NNODE), DEPTH, WL, WOW,SUMTZ,PRCNTER,
     * ITERA
  100 FORMAT ( 8F9.4, I4 )
      RETURN
      END
C
C
      SUBROUTINE PARAB (MXN,NNODE,X,TX,GW,RHS,PARA)
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXN), PARA(MXN),RHS(MXN)
      CONST = 0.5*GW/TX
      TL = X(NNODE) - X(1)
      C = (1./TL)*(RHS(NNODE)-CONST*TL*TL)
      DO I = 1 , NNODE
      PARA(I) = CONST*X(I)*X(I)+C*X(I)+RHS(1)
      END DO
      RETURN
      END
C
C
      SUBROUTINE CMPDEPTH (MXN,NNODE,X,RHS,DEPTH,WL )
C------- COMPUTES MAX. DEPTH AND LENGTH OF WIRE -------
C------- DEPTH = DEPTH OF WIRE FROM STRAIGHT LINE CONNECTING BOUNDARY.
C------- WL = LENGTH OF WIRE
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXN), RHS(MXN)
      WL = 0.
      SLOPE = (RHS(NNODE) - RHS(1))/(X(NNODE)-X(1))
      DEPTH = 0.
      DO I = 2 , NNODE
      Z0 = SLOPE*(X(I)-X(1)) + RHS(1)
      ZX = RHS(I) - Z0
      DEPTH = DMAX1 ( DEPTH, DABS(ZX) )
      DX = X(I) - X(I-1)
      DZ = RHS(I) - RHS(I-1)
      WL = WL + DSQRT ( DX*DX + DZ*DZ )
      END DO
      RETURN
      END
C
C
      SUBROUTINE TENSION (ND,MXE,MXN,NE,NODEX,X,RHS,TX,TS,SUMTZ)
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXN), RHS(MXN), TS(MXN), NODEX(MXE,ND)
C-------- TENSION IN EACH ELEMENT
      NNODE = NE + 1
      DO I = 1 , NNODE
      TS(I) = 0.
      END DO
      DO IEL = 1 , NE
      I = NODEX(IEL,1)
      J = NODEX(IEL,2)
      DX = X(J) - X(I)
      DZ = DABS (RHS(J) - RHS(I))
      DS = DSQRT( DX**2 + DZ**2 )
      TS(I) = TS(I) + TX*DS/DX
      TS(J) = TS(J) + TX*DS/DX
      END DO
      DO I = 2 , NE
      TS(I) = TS(I) / 2.
      END DO
C-------- UPSTREAM
      DX = X(2) - X(1)
      DZ = DABS ( RHS(2) - RHS(1) )
      SUMTZ = TX*DZ/DX
C-------- DWONSTREAM
      DX = X(NNODE) - X(NNODE-1)
      DZ = DABS ( RHS(NNODE) - RHS(NNODE-1) )
      SUMTZ = SUMTZ + TX*DZ/DX
      RETURN
      END
C
C
      SUBROUTINE MATRIX ( MXE,MXN,ND,NBW,NE,NNODE, GRAVITY,WDENSITY, 
     * TX,STIFF, NODEX, X, Z, A,RHS )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),X(MXN),Z(MXN),RHS(MXN),
     * STIFF(ND,ND),
     * A(MXN,NBW)  
      CONST = GRAVITY * WDENSITY / TX
      DO I = 1 , NNODE
      A(I,1) = 0.
      A(I,2) = 0.
      A(I,3) = 0.
      RHS(I) = 0.
      END DO
      DO IEL = 1 , NE
      I = NODEX(IEL,1)
      J = NODEX(IEL,2)
      DX = X(J) - X(I)
      DZ = Z(J) - Z(I)
      DSDX = DSQRT ( 1.D0 + (DZ/DX)**2 )
      FORCE = CONST * DSDX * DX / 2.
      EE = 1./DX
      STIFF(1,1) = -EE
      STIFF(1,2) =  EE
      STIFF(2,1) =  EE
      STIFF(2,2) = -EE
      A(I,2) = A(I,2) + STIFF(1,1)
      A(I,3) = A(I,3) + STIFF(1,2)
      A(J,1) = A(J,1) + STIFF(2,1)
      A(J,2) = A(J,2) + STIFF(2,2)
      RHS(I) = RHS(I) + FORCE
      RHS(J) = RHS(J) + FORCE
      END DO
      RETURN
      END
C
C
      SUBROUTINE INITIALZ ( MXN,NNODE,GW,TXMAX,X, Z1, ZN, Z )
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(MXN), Z(MXN)
      CONST = 0.5*GW/TXMAX
      TL = X(NNODE) - X(1)
      C = (1./TL)*(ZN-CONST*TL*TL)
      DO I = 1 , NNODE
      Z(I) = CONST*X(I)*X(I)+C*X(I)+Z1
      END DO
      RETURN
      END
C
C
      SUBROUTINE NEWZ ( MXN,NNODE,DOMAINL, RHS ,Z, DZMAXER )
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN), Z(MXN)
      DZMAXER = DABS ( RHS(1) - Z(1) ) / DOMAINL
      DO I = 1 , NNODE
      DZMAXER = DMAX1 ( DZMAXER, DABS(RHS(I)-Z(I))/DOMAINL  )
      Z(I) = (RHS(I) + Z(I))/2.
      END DO
      RETURN
      END
C
C
      SUBROUTINE INPUT (NE,DOMAINL,WDENSITY,TXMAX,TXMIN,TXSTEP, Z1,ZN)
      IMPLICIT REAL*8 ( A-H , O-Z )
      COMPLEX*16 UPSTRM, DWSTRM
      OPEN ( 1,FILE='WIRE.DAT', STATUS='OLD')
C-------- DOMAINL = LENGTH OF DOMAIN IN Z
C-------- WDENSITY = DENSITY OF WIRE AS GRAM PER CENTIMETER
C-------- TXMAX = POSSIBLE MAX. TENSION IN X-DIRECTION
C-------- YXMIN = POSSIBLE MIN. TENSION IN X-DIRECTION
C-------- TXSTEP = INCREMENTAL TENSION
C-------- Z1 = ELEVATION OF WIRE AT X=0
C-------- ZN = ELEVATION OF WIRE AT X=DOMAINL
C------- READING CONTROL PARAMETERS ------------------
      READ(1,*) NE
      READ(1,*) DOMAINL
      READ(1,*) WDENSITY
      READ(1,*) TXMAX,TXMIN,TXSTEP
      READ(1,*) Z1
      READ(1,*) ZN
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE ELEMENT ( DOMAINL,MXE,MXN,ND,NE,NNODE,NODEX,X  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),X(MXN)
C--------NE = NUMBER OF ELEMENT
      SLOPE = 0.
      A = 2.*(SLOPE-1.)
      B = 3.*(1.-SLOPE)
      DO I = 1 , NE
      NODEX(I,1) = I
      NODEX(I,2) = I+1
      END DO
      NNODE = NE + 1
      DX = DOMAINL / NE
      X(1) = 0.
      DO I = 2 , NNODE
      T = (I-1)*DX/DOMAINL
      X(I) = DOMAINL * T*(T*(A*T + B) + SLOPE)
      END DO
      RETURN
      END

C
C
      SUBROUTINE FORM ( MXN, NBW, NNODE, A, RHS, Z1, ZN )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION  A(MXN,NBW), RHS(MXN)
C-------UP STREAM SIDE
      A(1,2) = 1.
      A(1,3) = 0.
      RHS(1) = Z1
      RHS(2) = RHS(2) - Z1*A(2,1)
      A(2,1) = 0.
C-------DOWN STREAM SIDE
      A(NNODE,2) = 1.
      A(NNODE,1) = 0.
      RHS(NNODE) = ZN
      RHS(NNODE-1) = RHS(NNODE-1) - ZN*A(NNODE-1,3)
      A(NNODE-1,3) = 0.
      RETURN
      END
C
C
      SUBROUTINE SYSTEM ( MXN, NBW, NNODE, A , B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION  A(MXN,NBW) , B(MXN)
      B(1) = B(1) / A(1,2)
      A(1,2) = A(1,3) / A(1,2)
      DO I = 2 , NNODE
      P = A(I,2) - A(I,1) * A(I-1,2)
      A(I,2) = A(I,3) / P
      B(I) = ( B(I) - A(I,1)*B(I-1) ) / P
      END DO
C------ BACK SUBSTITUTION ----
      DO I = NNODE-1, 1,-1
      B(I) = B(I) - A(I,2) * B(I+1)
      END DO
      RETURN
      END
