sfellipi.red 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408
  1. module sfellipi; % Procedures and Rules for Elliptic Integrals.
  2. % Author: Lisa Temme, ZIB, October 1994
  3. algebraic <<
  4. %######################################################################
  5. %DESCENDING LANDEN TRANSFORMATION
  6. procedure landentrans(phi,alpha);
  7. begin scalar alpha_n!+1, alpha_n, phi_n!+1, phi_n, aNtoa0, pNtop0,
  8. a0toaN, p0topN;
  9. alpha_n := alpha;
  10. phi_n := phi;
  11. aNtoa0 := {alpha_n};
  12. pNtop0 := {phi_n};
  13. while alpha_n > 10^(-(Symbolic !:prec!:)) do
  14. <<
  15. alpha_n!+1:= asin(2/(1+cos(alpha_n)) -1);
  16. phi_n!+1 := phi_n + (atan(cos(alpha_n)*tan(phi_n)))
  17. + floor((floor(phi_n/(pi/2))+1)/2)*pi;
  18. aNtoa0 := alpha_n!+1.aNtoa0;
  19. pNtop0 := phi_n!+1.pNtop0;
  20. alpha_n := alpha_n!+1;
  21. phi_n := phi_n!+1
  22. >>;
  23. a0toaN := reverse(aNtoa0);
  24. p0topN := reverse(pNtop0);
  25. return list(p0topN, a0toaN)
  26. end;
  27. %######################################################################
  28. %VALUE OF EllipticF(phi,m)
  29. procedure F_function(phi,m);
  30. begin scalar alpha, bothlists, a0toaN, a1toaN, p0topN, phi_n, y,
  31. elptF;
  32. alpha := asin(sqrt(m));
  33. bothlists := landentrans(phi,alpha);
  34. a0toaN := PART(bothlists,2);
  35. a1toaN := REST(a0toaN);
  36. p0topN := PART(bothlists,1);
  37. phi_n := PART(reverse(p0topN),1);
  38. if phi = (pi/2)
  39. then
  40. elptF := K_function(m)
  41. else
  42. elptF :=
  43. phi_n *for each y in a1toaN PRODUCT(1/2)*(1+sin(y));
  44. return elptF
  45. end;
  46. %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  47. %EllipticF definition
  48. %====================
  49. operator EllipticF;
  50. EllipticFrules :=
  51. {
  52. EllipticF(~phi,0) => phi,
  53. EllipticF(i*~phi,0) => i*phi,
  54. EllipticF(~phi,1) => ln(sec(phi)+tan(phi)),
  55. EllipticF(i*~phi,1) => i*atan(sinh(phi)),
  56. EllipticF(~phi,~m) => Num_Elliptic(F_function,phi,m)
  57. when lisp !*rounded and numberp phi
  58. and numberp m
  59. };
  60. let EllipticFrules;
  61. %######################################################################
  62. %VALUE OF K(m)
  63. procedure K_function(m);
  64. begin scalar AGM, aN;
  65. AGM := AGM_function(1,sqrt(1-m),sqrt(m));
  66. aN := PART(AGM,2);
  67. return (pi / (2*aN));
  68. end;
  69. %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  70. %EllipticK definition
  71. %====================
  72. EllipticKrules :=
  73. {
  74. EllipticK(~m) => K_function(m) when lisp !*rounded
  75. and numberp m,
  76. EllipticK!'(~m) => K_function(1-m) when lisp !*rounded
  77. and numberp m
  78. };
  79. let EllipticKrules;
  80. %######################################################################
  81. %VALUE OF EllipticE(phi,m)
  82. procedure E_function(phi,m);
  83. begin scalar F, N, alpha, bothlists, a0toaN, p0topN, a1toaN, p1topN,
  84. sinalist, sinplist, b, s, blist, c, allz, w, z, allx,
  85. h, x, elptE;
  86. F := F_function(phi,m);
  87. alpha := asin(sqrt(m));
  88. bothlists := landentrans(phi,alpha);
  89. a0toaN := PART(bothlists, 2);
  90. p0topN := PART(bothlists, 1);
  91. a1toaN := REST(a0toaN);
  92. p1topN := REST(p0topN);
  93. N := LENGTH(a1toaN);
  94. sinalist := sin(a1toaN);
  95. sinplist := sin(p1topN);
  96. b := PART(sinalist,1);
  97. s := b;
  98. blist := for each c in rest sinalist collect << b := b*c >>;
  99. blist := s.blist;
  100. allz := 0;
  101. for w := 1:N do
  102. <<
  103. z := (1/(2^w))*PART(blist,w);
  104. allz := allz + z
  105. >>;
  106. allx := 0;
  107. for h := 1:N do
  108. <<
  109. x := (1/(2^h))*((PART(blist,h))^(1/2))
  110. * PART(sinplist,h);
  111. allx := allx + x
  112. >>;
  113. elptE := F * (1 - (1/2)*((sin(PART(a0toaN,1)))^2)*(1 + allz))
  114. + sin(PART(a0toaN,1))*allx ;
  115. return elptE;
  116. end;
  117. %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  118. %EllipticE(phi,m) definition
  119. %====================
  120. operator EllipticE;
  121. JacobiErules :=
  122. {
  123. EllipticE(0,~m) => 0,
  124. EllipticE(~phi,0) => phi,
  125. EllipticE(i*~phi,0) => i*phi,
  126. EllipticE(~phi,1) => sin(phi),
  127. EllipticE(i*~phi,1) => i*sinh phi,
  128. EllipticE(-~phi,~m) => -EllipticE(phi,m),
  129. EllipticE(~phi,-~m) => EllpiticE(phi,m),
  130. df(EllipticE(~phi,~m),~phi) => Jacobidn(phi,m)^2,
  131. df(EllipticE(~phi,~m),~m) =>
  132. m * (Jacobisn(phi,m) * Jacobicn(phi,m) * Jacobidn(phi,m)
  133. - EllipticE(phi,m) * Jacobicn(phi,m)^2) / (1-m^2)
  134. - m * phi * Jacobisn(phi,m)^2,
  135. EllipticE(~phi,~m) => Num_Elliptic(E_function,phi,m)
  136. when lisp !*rounded and numberp phi
  137. and numberp m,
  138. EllipticE(~m) => Num_Elliptic(E_function,pi/2,m)
  139. when lisp !*rounded and numberp m
  140. };
  141. let JacobiErules;
  142. %######################################################################
  143. %CALCULATING THE FOUR THETA FUNCTIONS
  144. %Theta 1 (often written H(u) - and has period 4K)
  145. %Theta 2 (often written H1(u) -and has period 4K)
  146. %Theta 3 (often written Theta1(u) - and has period 2K)
  147. %Theta 4 (often written Theta(u) - and has period 2K)
  148. procedure num_theta(a,u,m);
  149. begin scalar n, new, all, z, q, total;
  150. n := if a>2 then 1 else 0;
  151. new := 100; % To initiate loop
  152. all := 0;
  153. z := (pi*u)/(2*EllipticK(m));
  154. q := EXP(-pi*EllipticK(1-m)/EllipticK(m));
  155. while new > 10^(-(Symbolic !:prec!:)) do
  156. << new := if a =1 then
  157. ((-1)^n)*(q^(n*(n+1)))*sin((2*n+1)*z)
  158. else if a=2 then (q^(n*(n+1)))*cos((2*n+1)*z)
  159. else if a=3 then (q^(n*n))*cos(2*n*z)
  160. else if a=4 then ((-1)^n)*(q^(n*n))*cos(2*n*z);
  161. all := new + all;
  162. n := n+1
  163. >>;
  164. return if a > 2 then (1 + 2*all)
  165. else (2*(q^(1/4))*all);
  166. end;
  167. %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  168. %Theta Functions
  169. operator EllipticTheta;
  170. EllipticTHETArules :=
  171. {
  172. %Theta1rules
  173. %-----------
  174. EllipticTheta(1,~u,~m) =>
  175. Num_Elliptic(num_theta,1,u,m) when lisp !*rounded
  176. and numberp u and numberp m,
  177. EllipticTheta(1,-~u,~m) => -EllipticTheta(1,u,m),
  178. EllipticTheta(1,~u+EllipticK(~m),~m) => EllipticTheta(2,u,m),
  179. EllipticTheta(1,~u+(2*EllipticK(~m)),~m) =>
  180. -EllipticTheta(1,u,m),
  181. EllipticTheta(1,~u+i*EllipticK!'(~m),~m) =>
  182. i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  183. *EllipticTheta(4,u,m),
  184. EllipticTheta(1,~u+2*i*EllipticK!'(~m),~m) =>
  185. -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  186. *EllipticTheta(1,u,m),
  187. EllipticTheta(1,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
  188. (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  189. *EllipticTheta(3,u,m),
  190. EllipticTheta(1,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
  191. (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  192. *EllipticTheta(1,u,m),
  193. %Theta2rules
  194. %-----------
  195. EllipticTheta(2,~u,~m) =>
  196. Num_Elliptic(num_theta,2,u,m) when lisp !*rounded
  197. and numberp u and numberp m,
  198. EllipticTheta(2,-~u,~m) => EllipticTheta(2,u,m),
  199. EllipticTheta(2,~u+EllipticK(~m),~m) => -EllipticTheta(1,u,m),
  200. EllipticTheta(2,~u+(2*EllipticK(~m)),~m) =>
  201. -EllipticTheta(2,u,m),
  202. EllipticTheta(2,~u+i*EllipticK!'(~m),~m) =>
  203. (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  204. *EllipticTheta(3,u,m),
  205. EllipticTheta(2,~u+2*i*EllipticK!'(~m),~m) =>
  206. (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  207. *EllipticTheta(2,u,m),
  208. EllipticTheta(2,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
  209. -i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  210. *EllipticTheta(4,u,m),
  211. EllipticTheta(2,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
  212. -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  213. *EllipticTheta(2,u,m),
  214. %Theta3rules
  215. %-----------
  216. EllipticTheta(3,~u,~m) =>
  217. Num_Elliptic(num_theta,3,u,m) when lisp !*rounded
  218. and numberp u and numberp m,
  219. EllipticTheta(3,-~u,~m) => EllipticTheta(3,u,m),
  220. EllipticTheta(3,~u+EllipticK(~m),~m) => EllipticTheta(4,u,m),
  221. EllipticTheta(3,~u+(2*EllipticK(~m)),~m) =>
  222. EllipticTheta(3,u,m),
  223. EllipticTheta(3,~u+i*EllipticK!'(~m),~m) =>
  224. (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  225. *EllipticTheta(2,u,m),
  226. EllipticTheta(3,~u+2*i*EllipticK!'(~m),~m) =>
  227. (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  228. *EllipticTheta(3,u,m),
  229. EllipticTheta(3,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
  230. i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  231. *EllipticTheta(1,u,m),
  232. EllipticTheta(3,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
  233. (EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  234. *EllipticTheta(3,u,m),
  235. %Theta4rules
  236. %-----------
  237. EllipticTheta(4,~u,~m) =>
  238. Num_Elliptic(num_theta,4,u,m) when lisp !*rounded
  239. and numberp u and numberp m,
  240. EllipticTheta(4,-~u,~m) => EllipticTheta(4,u,m),
  241. EllipticTheta(4,~u+EllipticK(~m),~m) => EllipticTheta(3,u,m),
  242. EllipticTheta(4,~u+(2*EllipticK(~m)),~m)=>EllipticTheta(4,u,m),
  243. EllipticTheta(4,~u+i*EllipticK!'(~m),~m) =>
  244. i*(EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  245. *EllipticTheta(1,u,m),
  246. EllipticTheta(4,~u+2*i*EllipticK!'(~m),~m) =>
  247. -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  248. *EllipticTheta(4,u,m),
  249. EllipticTheta(4,~u+EllipticK(~m)+i*EllipticK!'(~m),~m) =>
  250. (EXP(-i*pi*0.5*u/EllipticK(m)))*(nome_q^(-1/2))
  251. *EllipticTheta(2,u,m),
  252. EllipticTheta(4,~u+2*EllipticK(~m)+2*i*EllipticK!'(~m),~m) =>
  253. -(EXP(-i*pi*u/EllipticK(m)))*(nome_q^-1)
  254. *EllipticTheta(4,u,m),
  255. %Error
  256. %-----
  257. EllipticTheta(~a,~u,~m) =>
  258. printerr ("In EllipticTheta(a,u,m); a = 1,2,3 or 4.")
  259. when numberp a
  260. and not(fixp a and a<5 and a>0)
  261. };
  262. let EllipticTHETArules;
  263. %######################################################################
  264. %CALCULATING ZETA
  265. procedure ZETA_function(u,m);
  266. begin scalar phi_list, clist, L, j, z, cn, phi_n;
  267. phi_list := PHI_function(1,sqrt(1-m),sqrt(m),u);
  268. clist := PART(AGM_function(1,sqrt(1-m),sqrt(m)),5);
  269. L := LENGTH(phi_list);
  270. j := 1;
  271. z := 0;
  272. while j < L do
  273. <<
  274. cn := PART(clist,L-j);
  275. phi_n := PART(phi_list,1+j);
  276. z := cn*sin(phi_n) + z;
  277. j := j+1
  278. >>;
  279. return z
  280. end;
  281. %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  282. %JacobiZETA definition
  283. %=====================
  284. operator JacobiZeta;
  285. JacobiZETArules :=
  286. {
  287. JacobiZeta(~u,0) => 0,
  288. JacobiZeta(~u,1) => tanh(u),
  289. JacobiZeta(-~u,~m) => -JacobiZeta(u,m),
  290. JacobiZeta(~u+~v,~m) => JacobiZeta(u,m) + JacobiZeta(v,m) -
  291. (m*Jacobisn(u,m)*Jacobisn(v,m)
  292. *Jacobisn(u+v,m)),
  293. JacobiZeta(~u+2*EllipticK(~m),m) => JacobiZeta(u,m),
  294. JacobiZeta(EllipticK(~m) - ~u,m) =>
  295. -JacobiZeta(EllipticK(m)+u,m),
  296. % JacobiZeta(~u,~m) => JacobiZeta(u - EllipticK(m),m) -
  297. % m * Jacobisn(u - EllipticK(m),m)
  298. % * Jacobicd(u - EllipticK(m),m),
  299. JacobiZeta(~u,~m) => Num_Elliptic(ZETA_function,u,m)
  300. when lisp !*rounded and numberp u
  301. and numberp m
  302. };
  303. let JacobiZETArules;
  304. %######################################################################
  305. >>;
  306. endmodule;
  307. end;