Changeset 5bd5df
 Timestamp:
 03/12/11 23:21:55 (4 years ago)
 Branches:
 master
 Children:
 f4a60f
 Parents:
 390a74
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

qdelliptic.lisp
r390a74 r5bd5df 696 696 ;; PI(n; phim) = integrate(1/sqrt(1m*sin(x)^2)/(1n*sin(x)^2), x, 0, phi) 697 697 ;; 698 ;; 699 ;; Carlson writes 700 ;; 701 ;; P(phi,k,n) = integrate((1+n*sin(t)^2)^(1)*(1k^2*sin(t)^2)^(1/2), t, 0, phi) 702 ;; = sin(phi)*Rf(cos(phi)^2, 1k^2*sin(phi)^2, 1) 703 ;;  n/3*sin(phi)^3*Rj(cos(phi)^2, 1k^2*sin(phi)^2, 1, 1+n*sin(phi)^2) 704 ;; 705 ;; Note that this definition as a different sign for the n parameter from A&S! 698 706 (defun ellipticpi (n phi m) 699 707 "Compute elliptic integral of the third kind: 700 708 701 709 PI(n; phim) = integrate(1/sqrt(1m*sin(x)^2)/(1n*sin(x)^2), x, 0, phi)" 702 ;; Note: Carlson's DRJ has n defined as the negative of the n given703 ;; in A&S.704 710 (let* ((precision (floatcontagion n phi m)) 705 711 (n (applycontagion n precision)) … … 709 715 (sinphi (sin phi)) 710 716 (cosphi (cos phi)) 711 (k (sqrt m)) 712 (k2sin (* ( 1 (* k sinphi)) 713 (+ 1 (* k sinphi))))) 714 ( (* sinphi (carlsonrf (expt cosphi 2) k2sin 1)) 717 (msin2 ( 1 (* m sinphi sinphi))) 718 ( (* sinphi (carlsonrf (expt cosphi 2) msin2 1)) 715 719 (* (/ nn 3) (expt sinphi 3) 716 (carlsonrj (expt cosphi 2) k2sin1720 (carlsonrj (expt cosphi 2) msin2 1 717 721 (+ 1 (* nn (expt sinphi 2)))))))) 
rttests.lisp
r390a74 r5bd5df 980 980 nil) 981 981 982 #+nil983 (rt:deftest oct.ellipticpi.19.6.2.d984 (loop for k from 0 below 100985 for n = (random 1d0)986 for epi = (ellipticpi ( n) (/ (floatpi n) 2) n)987 for true = (+ (/ (floatpi n) 4 (sqrt (+ 1 (sqrt n))))988 (/ (elliptick n) 2))989 for result = (checkaccuracy 53 epi true)990 when result991 append (list (list (list k n) result)))992 nil)993 994 995 #996 982 ;; ellipticpi(n, phi, 0) = 997 ;; atan h(sqrt(1n)*tan(phi))/sqrt(1n)n < 1983 ;; atan(sqrt(1n)*tan(phi))/sqrt(1n) n < 1 998 984 ;; atanh(sqrt(n1)*tan(phi))/sqrt(n1) n > 1 999 985 ;; tan(phi) n = 1 986 ;; 987 ;; These are easy to derive if you look at the integral: 988 ;; 989 ;; elliptipi(n, phi, 0) = integrate(1/(1n*sin(t)^2), t, 0, phi) 990 ;; 991 ;; and this can be easily integrated to give the above expressions for 992 ;; the different values of n. 1000 993 (rt:deftest oct.ellipticpi.n0.d 994 ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] 1001 995 (loop for k from 0 below 100 1002 996 for phi = (random (/ pi 2)) 1003 997 for n = (random 1d0) 1004 998 for epi = (ellipticpi n phi 0) 1005 for true = (/ (atan h(* (tan phi) (sqrt ( 1 n))))999 for true = (/ (atan (* (tan phi) (sqrt ( 1 n)))) 1006 1000 (sqrt ( 1 n))) 1007 for result = (checkaccuracy 5 3epi true)1001 for result = (checkaccuracy 50 epi true) 1008 1002 unless (eq nil result) 1009 1003 append (list (list (list k n phi) result))) … … 1013 1007 (loop for k from 0 below 100 1014 1008 for phi = (random (/ pi 2)) 1015 for epi = (ellipticpi 0phi 0)1009 for epi = (ellipticpi 1 phi 0) 1016 1010 for true = (tan phi) 1017 for result = (checkaccuracy 53 epi true)1011 for result = (checkaccuracy 43 epi true) 1018 1012 unless (eq nil result) 1019 1013 append (list (list (list k phi) result))) … … 1027 1021 for true = (/ (atanh (* (tan phi) (sqrt ( n 1)))) 1028 1022 (sqrt ( n 1))) 1029 for result = (checkaccuracy 52epi true)1023 for result = (checkaccuracy 49 epi true) 1030 1024 ;; Not sure if this formula holds when atanh gives a complex 1031 1025 ;; result. Wolfram doesn't say … … 1035 1029 1036 1030 (rt:deftest oct.ellipticpi.n0.q 1031 ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] 1037 1032 (loop for k from 0 below 100 1038 1033 for phi = (random (/ +pi+ 2)) 1039 1034 for n = (random #q1) 1040 1035 for epi = (ellipticpi n phi 0) 1041 for true = (/ (atan h(* (tan phi) (sqrt ( 1 n))))1036 for true = (/ (atan (* (tan phi) (sqrt ( 1 n)))) 1042 1037 (sqrt ( 1 n))) 1043 for result = (checkaccuracy 2 12epi true)1038 for result = (checkaccuracy 208 epi true) 1044 1039 unless (eq nil result) 1045 1040 append (list (list (list k n phi) result))) … … 1049 1044 (loop for k from 0 below 100 1050 1045 for phi = (random (/ +pi+ 2)) 1051 for epi = (ellipticpi 0phi 0)1046 for epi = (ellipticpi 1 phi 0) 1052 1047 for true = (tan phi) 1053 for result = (checkaccuracy 2 12epi true)1048 for result = (checkaccuracy 205 epi true) 1054 1049 unless (eq nil result) 1055 1050 append (list (list (list k phi) result))) … … 1063 1058 for true = (/ (atanh (* (tan phi) (sqrt ( n 1)))) 1064 1059 (sqrt ( n 1))) 1065 for result = (checkaccuracy 20 9epi true)1060 for result = (checkaccuracy 208 epi true) 1066 1061 ;; Not sure if this formula holds when atanh gives a complex 1067 1062 ;; result. Wolfram doesn't say … … 1069 1064 append (list (list (list k n phi) result))) 1070 1065 nil) 1071 #
Note: See TracChangeset
for help on using the changeset viewer.