Sources Esope | Cast3M

cylwrk
C CYLWRK    SOURCE    GF238795  18/02/01    21:15:05     9724           
      SUBROUTINE CYLWRK(RF,ZF,KFLAG,U,UR,UZ,URR,URZ,UZZ,ITYPE,
     $NPN,EDICON,ICOO,ITRAVA)
C
C      EVALUATE POTENTIAL AND ITS DERIVATIVES FOR CYLINDER COORDINATES
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C      POTENTIAL VALUES
      SEGMENT EDICON
        INTEGER  KSTRT, KSTEP, NMIR, IS
        REAL*8  CROT, SROT, SYMFCT
        LOGICAL LREAL, LIMAG
      ENDSEGMENT
      SEGMENT ICOO
        REAL*8 X(MV),Y(MV),P(MV),WNODE(MV)
        INTEGER LISVO(MV)
      ENDSEGMENT
      SEGMENT ITRAVA
        REAL*8 KENN(M42,2),SIGMA(M42),DELRHO(M42),C(M50,M50)
        REAL*8 AK(M50),UM(M50),RM(M50)
       INTEGER IL(M50)
      ENDSEGMENT
C       COMMON // NEU, NPN, NPU
C
C      CONSTANTS DESCRIBING SYMMETRY CONDITIONS
C
C
C      SOURCE TERMS
C
C      TABLES AND WORKING STORE FOR POTENTIAL AND FIELD EDIT
C
C-----------------------------------------------------------------------
C
      M50=RM(/1)
      NPU = MIN(NPN,19)
      KEND = KSTRT + KSTEP * (NPU - IS - 2)
      IFLAG = 0
      J = 0
C
      DO 90 JJ = 1, NPN
C
      IF(WNODE(JJ) .EQ. 0.0) GO TO 90
      RN  = X(JJ)
      ZN  = Y(JJ)
      UN  = P(JJ)
C
C      SET UP HARMONIC FUNCTIONS FOR LOCATION (RN,ZN)
C      NOTE.         U(SUBSCRIPT M) IS STORED IN UM(M+3)
C                    R(SUBSCRIPT M) IS STORED IN RM(M+3)
C                    THE DO LOOP INDEX K IS (M+1)
C
      UM(3) = 1.0
      IF(IS .NE. 0) UM(3) = RN
      RM(3) = FLOAT(IS+IS) * ZN
C
      DO 40 K = 2, KEND
      UM(K+2) = (RN*RM(K+1) + ZN*UM(K+1))
      RM(K+2) = (ZN*RM(K+1) - RN*UM(K+1)) * FLOAT(K+IS+IS-1) / FLOAT(K)
   40 CONTINUE
C
C      SET UP COEFFICIENT MATRIX FOR FIT
C
      J = J + 1
      I = IS
      C(J,1) = WNODE(JJ)
C
      IF(KEND .LT. KSTRT) GO TO 80
C
      DO 70 I = IS+1, IS+1+(KEND-KSTRT)/KSTEP
      K = KSTRT + (I-IS-1)*KSTEP
c      DO 70 K = KSTRT, KEND, KSTEP
c      I = I + 1
      C(J,I) = WNODE(JJ) * UM(3) * UM(K+2)
   70 CONTINUE
C
   80 C(J,I+1) = WNODE(JJ) * UN
C
   90 CONTINUE
C      WRITE(6,FMT= ' ('' I ET J MAX '',/,(2I6))') I,J
C      WRITE(6,FMT= ' ('' KSTRT KSTEP KEND ITYPE  '',/,(4I5))') KSTRT,
C     $KSTEP,KEND,ITYPE
C
C-----------------------------------------------------------------------
C      PERFORM LEAST-SQUARES FIT
C
      CALL GLSQ(C,AK,IL,J,I,ALFA,1.0D-10,1.0D-10,M50)
C
C-----------------------------------------------------------------------
C      COMPUTE POTENTIAL AND DERIVATIVES
C
      UM(1) = 0.0
      UM(2) = 0.0
      UM(3) = 1.0
      IF(IS .NE. 0) UM(3) = RF
      RM(1) = 0.0
      RM(2) = FLOAT(IS)
      RM(3) = 2.0 * RM(2) * ZF
C
      DO 110 K = 2, KEND
      UM(K+2) = (RF*RM(K+1) + ZF*UM(K+1))
      RM(K+2) = (ZF*RM(K+1) - RF*UM(K+1)) * FLOAT(K+IS+IS-1) / FLOAT(K)
  110 CONTINUE
C
      IFLAG = KFLAG
      U   = AK(1) * FLOAT(IS)
      UR  = 0.0
      UZ  = 0.0
      URZ = 0.0
      UZZ = 0.0
      I   = IS
C
      IF(KEND .LT. KSTRT) GO TO 150
C
      DO 120 K = KSTRT, KEND, KSTEP
      I = I + 1
      U   = U   + AK(I) * UM(3) * UM(K+2)
      XK  = AK(I) * FLOAT(K+IS+IS-1)
      UR  = UR  + XK * RM(K+1)
      UZ  = UZ  + XK * UM(K+1)
      XK  = XK * FLOAT(K+IS+IS-2)
      URZ = URZ + XK * RM(K)
      UZZ = UZZ + XK * UM(K)
  120 CONTINUE
C
  150 KFLAG = IFLAG
      URR = - UZZ
      RETURN
C***********************************************************************
      END