github.com/gopherd/gonum@v0.0.4/mathext/ell_carlson_test.go (about)

     1  // Copyright ©2017 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package mathext
     6  
     7  import (
     8  	"math"
     9  	"testing"
    10  
    11  	"math/rand"
    12  )
    13  
    14  // Testing EllipticF (and EllipticRF) using the addition theorems from http://dlmf.nist.gov/19.11.i
    15  func TestEllipticF(t *testing.T) {
    16  	t.Parallel()
    17  	const tol = 1.0e-14
    18  	rnd := rand.New(rand.NewSource(1))
    19  
    20  	// The following EllipticF(pi/3,m), m=0.1(0.1)0.9 was computed in Maxima 5.38.0 using Bigfloat arithmetic.
    21  	vF := [...]float64{
    22  		1.0631390181954904767742338285104637431858016483079,
    23  		1.0803778062523490005579242592072579594037132891908,
    24  		1.0991352230920430074586978843452269008747645822123,
    25  		1.1196949183404746257742176145632376703505764745654,
    26  		1.1424290580457772555013955266260457822322036529624,
    27  		1.1678400583161860445148860686430780757517286094732,
    28  		1.1966306515644649360767197589467723191317720122309,
    29  		1.2298294422249382706933871574135731278765534034979,
    30  		1.2690359140762658660446752406901433173504503955036,
    31  	}
    32  	phi := math.Pi / 3
    33  	for m := 1; m <= 9; m++ {
    34  		mf := float64(m) / 10
    35  		delta := math.Abs(EllipticF(phi, mf) - vF[m-1])
    36  		if delta > tol {
    37  			t.Fatalf("EllipticF(pi/3,m) test fail for m=%v", mf)
    38  		}
    39  	}
    40  
    41  	for test := 0; test < 100; test++ {
    42  		alpha := rnd.Float64() * math.Pi / 4
    43  		beta := rnd.Float64() * math.Pi / 4
    44  		for mi := 0; mi < 9999; mi++ {
    45  			m := float64(mi) / 10000
    46  			Fa := EllipticF(alpha, m)
    47  			Fb := EllipticF(beta, m)
    48  			sina, cosa := math.Sincos(alpha)
    49  			sinb, cosb := math.Sincos(beta)
    50  			tan := (sina*math.Sqrt(1-m*sinb*sinb) + sinb*math.Sqrt(1-m*sina*sina)) / (cosa + cosb)
    51  			gamma := 2 * math.Atan(tan)
    52  			Fg := EllipticF(gamma, m)
    53  			delta := math.Abs(Fa + Fb - Fg)
    54  			if delta > tol {
    55  				t.Fatalf("EllipticF test fail for m=%v, alpha=%v, beta=%v", m, alpha, beta)
    56  			}
    57  		}
    58  	}
    59  }
    60  
    61  // Testing EllipticE (and EllipticRF, EllipticRD) using the addition theorems from http://dlmf.nist.gov/19.11.i
    62  func TestEllipticE(t *testing.T) {
    63  	t.Parallel()
    64  	const tol = 1.0e-14
    65  	rnd := rand.New(rand.NewSource(1))
    66  
    67  	// The following EllipticE(pi/3,m), m=0.1(0.1)0.9 was computed in Maxima 5.38.0 using Bigfloat arithmetic.
    68  	vE := [...]float64{
    69  		1.0316510822817691068014397636905610074934300946730,
    70  		1.0156973658341766636288643556414001451527597364432,
    71  		9.9929636467826398814855428365155224243586391115108e-1,
    72  		9.8240033979859736941287149003648737502960015189033e-1,
    73  		9.6495145764299257550956863602992167490195750321518e-1,
    74  		9.4687829659158090935158610908054896203271861698355e-1,
    75  		9.2809053417715769009517654522979827392794124845027e-1,
    76  		9.0847044378047233264777277954768245721857017157916e-1,
    77  		8.8785835036531301307661603341327881634688308777383e-1,
    78  	}
    79  	phi := math.Pi / 3
    80  	for m := 1; m <= 9; m++ {
    81  		mf := float64(m) / 10
    82  		delta := math.Abs(EllipticE(phi, mf) - vE[m-1])
    83  		if delta > tol {
    84  			t.Fatalf("EllipticE(pi/3,m) test fail for m=%v", mf)
    85  		}
    86  	}
    87  
    88  	for test := 0; test < 100; test++ {
    89  		alpha := rnd.Float64() * math.Pi / 4
    90  		beta := rnd.Float64() * math.Pi / 4
    91  		for mi := 0; mi < 9999; mi++ {
    92  			m := float64(mi) / 10000
    93  			Ea := EllipticE(alpha, m)
    94  			Eb := EllipticE(beta, m)
    95  			sina, cosa := math.Sincos(alpha)
    96  			sinb, cosb := math.Sincos(beta)
    97  			tan := (sina*math.Sqrt(1-m*sinb*sinb) + sinb*math.Sqrt(1-m*sina*sina)) / (cosa + cosb)
    98  			gamma := 2 * math.Atan(tan)
    99  			Eg := EllipticE(gamma, m)
   100  			delta := math.Abs(Ea + Eb - Eg - m*sina*sinb*math.Sin(gamma))
   101  			if delta > tol {
   102  				t.Fatalf("EllipticE test fail for m=%v, alpha=%v, beta=%v", m, alpha, beta)
   103  			}
   104  		}
   105  	}
   106  }