elipt.f90

   1 MODULE ELIPT_MOD
   2
   3   IMPLICIT NONE
   4
   5 CONTAINS
   6
   7 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   8 !                Heuman's function  (KP!!2 = 1-K!!2)                   !
   9 !                                                                      !
  10 !                     (PI/2)* LAMBDA0(PSI,K) =                         !
  11 !                                                                      !
  12 !          = E(K)*F(PSI,KP) + K(K)*( E(PSI,KP) - F(PSI,KP) )           !
  13 !                                                                      !
  14 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  15
  16       DOUBLE PRECISION FUNCTION LAMBDA0(PSI,K)
  17
  18 !!!!! variables 'dummy' y variables internas !!!!!!!!!!!!!!!!!!!!!!!!!!!
  19
  20       DOUBLE PRECISION &
  21                    PSI,K,KFUN,EFUN,KP,FINTEL,EINTEL,PI,DPI
  22
  23 !!!!! COMMON variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  24
  25
  26       PI  = DATAN(1.D0)*4.D0
  27       DPI = 2.D0*PI
  28
  29       CALL ELIPT0(K,KFUN,EFUN)
  30
  31       KP = DSQRT(1.D0-K**2)
  32
  33       CALL ELIPT(PSI,KP,FINTEL,EINTEL)
  34
  35       LAMBDA0 = (2.D0/PI)*( EFUN*FINTEL + KFUN*( EINTEL-FINTEL ) )
  36
  37       RETURN
  38       END FUNCTION LAMBDA0
  39
  40 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  41 !     Incomplete elliptic integrals
  42 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  43
  44       SUBROUTINE ELIPT(PSI,K,FINTEL,EINTEL)
  45
  46 !!!!! variables 'dummy' y variables internas !!!!!!!!!!!!!!!!!!!!!!!!!!!
  47
  48       DOUBLE PRECISION  PSI,K,FINTEL,EINTEL,S,C
  49       INTEGER      IERR
  50
  51 !!!!! function names !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  52
  53       DOUBLE PRECISION DRF,DRD
  54
  55 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  56
  57       S = DSIN(PSI)
  58       C = DCOS(PSI)
  59
  60       FINTEL = S*DRF(C**2,1.D0-(K*S)**2,1.D0,IERR)
  61
  62       IF (IERR.NE.0) THEN
  63          WRITE(6,*) 'STOP EN FINTEL','IERR =',IERR
  64          STOP
  65       END IF
  66
  67       EINTEL = FINTEL  &
  68              - (K**2/3.D0)*S**3*DRD(C**2,1.D0-(K*S)**2,1.D0,IERR)
  69
  70       IF (IERR.NE.0) THEN
  71          WRITE(6,*) 'STOP EN EINTEL','IERR =',IERR
  72          STOP
  73       END IF
  74
  75       RETURN
  76       END SUBROUTINE ELIPT
  77
  78 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  79 !     Complete elliptic integrals
  80 !     poner los limites en K = 0
  81 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  82
  83       SUBROUTINE ELIPT0(K,KFUN,EFUN)
  84
  85 !!!!! variables 'dummy' y variables internas !!!!!!!!!!!!!!!!!!!!!!!!!!!
  86
  87       DOUBLE PRECISION K,KFUN,EFUN
  88       INTEGER      IERR
  89
  90 !!!!! function names !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  91
  92       DOUBLE PRECISION DRF,DRD
  93
  94 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  95
  96       KFUN = DRF(0.D0,1.D0-K**2,1.D0,IERR)
  97
  98       IF (IERR.NE.0) THEN
  99          WRITE(6,*) 'STOP IN KFUN','IERR =',IERR
 100          STOP
 101       END IF
 102
 103       EFUN = KFUN  &
 104            - (K**2/3.D0)*DRD(0.D0,1.D0-K**2,1.D0,IERR)
 105
 106       IF (IERR.NE.0) THEN
 107          WRITE(6,*) 'STOP IN EFUN','IERR =',IERR
 108          STOP
 109       END IF
 110
 111       RETURN
 112       END SUBROUTINE ELIPT0
 113
 114       SUBROUTINE ELIPTK0(K,KFUN)
 115
 116 !!!!! variables 'dummy' y variables internas !!!!!!!!!!!!!!!!!!!!!!!!!!!
 117
 118       DOUBLE PRECISION K,KFUN
 119       INTEGER      IERR
 120
 121 !!!!! function names !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 122
 123       DOUBLE PRECISION DRF
 124
 125 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 126
 127       KFUN = DRF(0.D0,1.D0-K**2,1.D0,IERR)
 128
 129       IF (IERR.NE.0) THEN
 130          WRITE(6,*) 'STOP IN KFUN','IERR =',IERR
 131          STOP
 132       END IF
 133
 134       RETURN
 135       END SUBROUTINE ELIPTK0
 136
 137 END MODULE ELIPT_MOD