Numerical Computation of the Elliptic Functions 

by Michael W. Pashea



Listing One



/***********************************************************************

  Subroutines for numerical calculation of the elliptic functions sn(u,k),

  cn(u,k), and dn(u,k). Also included are subroutines to calculate the 

  first and second types of elliptic integrals, F(phi,k), K(k) and E(phi,k).

  Michael W. Pashea 11-1-04  

***********************************************************************/

#include <stdio.h>

#include <math.h>



// Define PI, Tolerance and Maximum Iterations

#define PI 3.14159265358979

#define TOL 0.0000000001

#define MAX_ITERATIONS 10



// Function prototypes - normally this would be in a header file

double find_m( double angle_in_rad );

double k_next( double k );

double phi_next( double phi, double k );

double sn( double u, double k );

double cn( double u, double k );

double dn( double u, double k );

double F_incomplete( double phi_in_rad, double k);

double K_complete( double k);

double E_incomplete( double phi_in_rad, double k);

double E_complete( double k);



// Main program - Just print out the data from the examples

void main(){

    printf("The elliptic sine sn(1.235,0.5) ");

    printf("evaluates to: %1.13lf \n", sn(1.235, 0.5) );

    printf("The elliptic cosine cn(1.235,0.5) ");

    printf("evaluates to: %1.13lf \n", cn(1.235, 0.5) );

    printf("The elliptic delta dn(1.235,0.5) ");

    printf("evaluates to: %1.13lf \n\n", dn(1.235, 0.5) );

    printf("The incomplete Elliptic integral F(PI/3,0.35) ");

    printf("evaluates to: %1.13lf \n", F_incomplete(PI/3, 0.35) );

    printf("The complete Elliptic integral K(0.5) ");

    printf("evaluates to: %1.13lf \n\n", K_complete(0.5) );

    printf("The complete Elliptic integral E(0.3) ");

    printf("evaluates to: %1.13lf \n", E_complete(0.3) );

    printf("The incomplete Elliptic integral E(PI/6, 0.5) ");

    printf("evaluates to: %1.13lf \n", E_incomplete(PI/6, 0.5) );

}

// 

double find_m( double angle_in_rad ){

    return( floor( angle_in_rad/PI + 0.5));

}   

double k_next( double k ){

    k = ( 1 - sqrt( 1-k*k ))/( 1 + sqrt( 1-k*k ));

    return( k );

}

double phi_next( double phi, double k ){

    double m;

    m = find_m(phi); 

    k = k_next( k );

    phi = phi + (m * PI) + atan( ((1-k)/(1+k))*tan( phi-(m * PI) ) ); 

    return( phi );

}

// sn(u,k), cn(u,k) and dn(u,k) as demonstrated in Example 1

// These could be easily combined into one subroutine 

double sn( double u, double k ){

    int i,Nmax;

    double kvalue[MAX_ITERATIONS];

    Nmax = 1;

    for(i=0; i<MAX_ITERATIONS; i++){

        kvalue[i] = 0.0;

    }

    kvalue[0] = k;

    while( k > TOL ){

        k = k_next(k);

        kvalue[Nmax] = k;

        u *= 2/(1+k);

        Nmax++;

    }

    for (i=Nmax-1; i>0; i--){

        u = (u + asin( kvalue[i] * sin(u) ))/2;

    }

    return(sin(u));

}

double cn( double u, double k ){

    int i,Nmax;

    double kvalue[MAX_ITERATIONS];

    Nmax = 1;

    for(i=0; i<MAX_ITERATIONS; i++){

        kvalue[i] = 0.0;

    }

    kvalue[0] = k;

    while( k > TOL ){

        k = k_next(k);

        kvalue[Nmax] = k;

        u *= 2/(1+k);

        Nmax++;

    }

    for (i=Nmax-1; i>0; i--){

        u = (u + asin( kvalue[i] * sin(u) ))/2;

    }

    return(cos(u));

}

double dn( double u, double k ){

    int i,Nmax;

    double kvalue[MAX_ITERATIONS];

    Nmax = 1;

    for(i=0; i<MAX_ITERATIONS; i++){

        kvalue[i] = 0.0;

    }

    kvalue[0] = k;

    while( k > TOL ){

        k = k_next(k);

        kvalue[Nmax] = k;

        u *= 2/(1+k);

        Nmax++;

    }

    for (i=Nmax-1; i>0; i--){

        u = (u + asin( kvalue[i] * sin(u) ))/2;

    }

    return(sqrt(1-(kvalue[0]*kvalue[0]*sin(u))));

}

// Incomplete elliptic integral of the first kind as in Example 4

double F_incomplete( double phi_in_rad, double k){

    double F;

    F = 1.0;

    while ( k > TOL ){

        phi_in_rad = phi_next(phi_in_rad, k);

        k = k_next( k );

        F *= (1 + k)/2;

    }

   return (F * phi_in_rad);

}

// Complete elliptic integral of the first kind as in Example 2

double K_complete( double k){

    double K;

    K = PI/2;

    while ( k > TOL ){

        k = k_next( k );

        K *= (1 + k);

    }

    return (K);

}

// Incomplete elliptic integral of the second kind. Discussed in 

// Example 3, but not calculated. 

double E_incomplete( double phi_in_rad, double k ){

    double E, F;

    double a, aterm, asum;

    double bterm, bsum;



    E = 1.0;

    F = F_incomplete(phi_in_rad, k );



    a = (k*k)/2;

    aterm = 1.0;

    asum = 1.0;

    bterm = k;

    bsum = 0.0;



    while ( k > TOL ){

        bterm = bterm/k;

        phi_in_rad = phi_next(phi_in_rad, k);

        k = k_next( k );

        aterm *= (aterm*k)/2;

        asum += aterm;              

        bterm *= k/(1+k);

        bsum+= bterm*sin(phi_in_rad);

    }

    E = (1-a*asum)*F+bsum;

    return (E);

}

// Complete elliptic integral of the second kind as in Example 3

double E_complete( double k ){

    double E;

    double a, aterm, asum;



    E = PI/2;

    a = (k*k)/2;

    aterm = 1.0;

    asum = 1.0;



    while ( k > TOL ){

        k = k_next( k );

        E *= (1 + k);

        aterm *= (aterm*k)/2;

        asum += aterm;

    }

    E = (1-a*asum)*E;

    return (E);

}







	   



4



