C Aproximace eliptickych integralu.
C Podle: Abramowicz,Stegun: Handbook of matehematical functions, 1964
C



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C Uplny elipticky integral prvniho druhu.
C
C        / 1
C K(m) = |    [(1 - t**2)*(1 - m*t**2)]**(-1/2) dt
C      0 /
C
C        / pi/2
C K(m) = |    (1 - m*sin(t)**2)**(-1/2) dt
C      0 /
C
C Presnost je pro 0 <= M < 1.0 je vetsi jak 2E-8

      REAL FUNCTION K1(M)
      REAL M
 
      INTEGER I
      DOUBLE PRECISION M1,A(0:4),B(0:4),SA,SB

      DATA A/ 1.38629 436112,
     .        0.09666 344259,
     .        0.03590 092383,
     .        0.03742 563713,
     .        0.01451 196212 /
      DATA B/ 0.50000 000000,
     .        0.12498 593597,
     .        0.06880 248576,
     .        0.03328 355346,
     .        0.00441 787012 /

      M1 = 1.0 - M
      SA = A(4)
      DO I = 3,0,-1
        SA = M1*SA + A(I)
      ENDDO
      SB = B(4)
      DO I = 3,0,-1
        SB = M1*SB + B(I)
      ENDDO

      K1 = SA + SB*LOG(1.0/M1)

      END


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C Uplny elipticky integral druheho druhu.
C
C        / 1
C K(m) = |    [(1 - t**2)*(1 - m*t**2)]**(1/2) dt
C      0 /
C
C        / pi/2
C K(m) = |    (1 - m*sin(t)**2)**(1/2) dt
C      0 /
C
C Presnost je pro 0 <= M < 1.0 je vetsi jak 2E-8

      REAL FUNCTION E(M)
      REAL M
 
      INTEGER I
      DOUBLE PRECISION M1,A(0:4),B(0:4),SA,SB

      DATA A/ 1.00000 000000,
     .        0.44325 141463,
     .        0.06260 601220,
     .        0.04757 383546,
     .        0.01736 506451 /
      DATA B/ 0.00000 000000,
     .        0.24998 368310,
     .        0.09200 180037,
     .        0.04069 697526,
     .        0.00526 449639 /

      M1 = 1.0 - M
      SA = A(4)
      DO I = 3,0,-1
        SA = M1*SA + A(I)
      ENDDO
      SB = B(4)
      DO I = 3,0,-1
        SB = M1*SB + B(I)
      ENDDO

      E = SA + SB*LOG(1.0/M1)

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL FUNCTION ELIPT1( PHI,ALPHA )
C Vypocet eliptickeho integralu prvniho druhu
C                       / phi
C   F( phi \ aplha ) = |   (1 - sin^2 alpha * sin^2 theta)^(-1/2) dtheta 
C                    0 /
C pomoci "procesu aritmeticko-geometrickeho prumeru" popsaneho
C v Abramowitz,Stegung: Handbook of special functions, 1965
C
C Vstupem jsou uhly v radianech.

      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415 92653 58979 32384 62643 )

      REAL PHI,ALPHA
      DOUBLE PRECISION A,B,A1,B1,X,EPS,F1,F2
      INTEGER N
      DATA EPS/1E-8/

      A = 1.0
      B = COS(ALPHA)
      N = 0
      F1 = PHI

10    IF( ABS(A-B).LT.EPS .OR. N.GT.100 ) GOTO 20
        N = N + 1     
        X = B/A
        F2 = X*TAN(F1)
        F2 = ATAN(F2)
        IF( F2.LT.0.0 ) F2 = F2 + PI 
        F1 = F2 + F1 + AINT(F1/PI)*PI
        A1 = A
        B1 = B
        A = 0.5*(A1 + B1)
        B = SQRT(A1*B1)
c        print *,N,A,B,F1/2**N/A
      GOTO 10
20    N = 2**N
      ELIPT1 = F1/N/A
      END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL FUNCTION ELIPT2( PHI,ALPHA )
C POZOR!!!!   Nefunguje!
C Vypocet eliptickeho integralu druheho druhu
C                       / phi
C   F( phi \ aplha ) = |   (1 - sin^2 alpha * sin^2 theta)^(1/2) dtheta 
C                    0 /
C pomoci "procesu aritmeticko-geometrickeho prumeru" popsaneho
C v Abramowitz,Stegung: Handbook of special functions, 1965
C
C Vstupem jsou uhly v radianech.

      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415 92653 58979 32384 62643 )

      REAL PHI,ALPHA,K1,E,M
      DOUBLE PRECISION A,B,C,A1,B1,N2,X,EPS,F1,F2,FF,EK,EE
      DOUBLE PRECISION SC2,SE
      INTEGER N
      EXTERNAL E,K1
      DATA EPS/1E-8/

      M = SIN(ALPHA)**2
      print *,M,E(M),K1(M)
      EK = E(M)/K1(M)
      A = 1.0
      B = COS(ALPHA)
      C = SIN(ALPHA)
      SC2 = C**2
      SE = 0.0
      N = 0
      N2 = 1.0
      F1 = PHI

10    IF( ABS(A-B).LT.EPS .OR. N.GT.100 ) GOTO 20
        N = N + 1  
        N2 = 2.0*N2   
        X = B/A
C vypocet elipt. integralu 1 druhu
        F2 = X*TAN(F1)
        F2 = ATAN(F2)
        IF( F2.LT.0.0 ) F2 = F2 + PI 
        F1 = F2 + F1 + AINT(F1/PI)*PI
        FF = F1/N2/A
C vypocet uplneho elipt. integralu 1 druhu
c        K = PI/2.0/A
C vypocet uplneho elipt. integralu 2 druhu
c        E = K*(1.0 - 0.5*SC2)
c        print *,E,SC2
C vypocet elipt. integralu 2 druhu
        EE = SE + EK*FF
        A1 = A
        B1 = B
        A = 0.5*(A1 + B1)
        B = SQRT(A1*B1)
        C = -0.5*(A - B)
        SC2 = SC2 + N2*C**2
        SE = SE + C*SIN(F1)
        print *,N,real(C),real(EK),real(FF),real(EE)
      GOTO 10
20    ELIPT2 = EE

      END