Changeset e2d8d6
 Timestamp:
 03/12/11 03:57:40 (4 years ago)
 Branches:
 master
 Children:
 147fa2
 Parents:
 5c21f1
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

qdelliptic.lisp
r5f6c62 re2d8d6 73 73 (doublefloat 'doublefloat) 74 74 (qdreal 'qdreal)))) 75 (if complexp 76 (if (eq maxtype 'qdreal) 77 'qdcomplex 78 `(cl:complex ,maxtype)) 79 maxtype))) 80 75 maxtype)) 76 77 (defun applycontagion (number precision) 78 (etypecase number 79 ((or cl:real qdreal) 80 (coerce number precision)) 81 ((or cl:complex qdcomplex) 82 (complex (coerce (realpart number) precision) 83 (coerce (imagpart number) precision))))) 84 81 85 ;;; Jacobian elliptic functions 82 86 … … 284 288 285 289 Rf(x, y, z) = 1/2*integrate((t+x)^(1/2)*(t+y)^(1/2)*(t+z)^(1/2), t, 0, inf)" 286 (let* ((xn x) 287 (yn y) 288 (zn z) 290 (let* ((precision (floatcontagion x y z)) 291 (xn (applycontagion x precision)) 292 (yn (applycontagion y precision)) 293 (zn (applycontagion z precision)) 289 294 (a (/ (+ xn yn zn) 3)) 290 295 (epslon (/ (max (abs ( a xn)) … … 335 340 336 341 Rd(x,y,z) = integrate(3/2*(t+x)^(1/2)*(t+y)^(1/2)*(t+z)^(3/2), t, 0, inf)" 337 (let* ((xn x) 338 (yn y) 339 (zn z) 342 (let* ((precision (floatcontagion x y z)) 343 (xn (applycontagion x precision)) 344 (yn (applycontagion y precision)) 345 (zn (applycontagion z precision)) 340 346 (a (/ (+ xn yn (* 3 zn)) 5)) 341 347 (epslon (/ (max (abs ( a xn)) … … 350 356 (loop while (> (* power4 epslon) (abs an)) 351 357 do 352 353 354 355 356 357 358 359 360 361 362 363 364 365 358 (setf xnroot (sqrt xn)) 359 (setf ynroot (sqrt yn)) 360 (setf znroot (sqrt zn)) 361 (setf lam (+ (* xnroot ynroot) 362 (* xnroot znroot) 363 (* ynroot znroot))) 364 (setf sigma (+ sigma (/ power4 365 (* znroot (+ zn lam))))) 366 (setf power4 (* power4 1/4)) 367 (setf xn (* (+ xn lam) 1/4)) 368 (setf yn (* (+ yn lam) 1/4)) 369 (setf zn (* (+ zn lam) 1/4)) 370 (setf an (* (+ an lam) 1/4)) 371 (incf n)) 366 372 ;; c1=3/14,c2=1/6,c3=9/88,c4=9/22,c5=3/22,c6=9/52,c7=3/26 367 373 (let* ((xndev (/ (* ( a x) power4) an)) … … 386 392 (* 45/272 ee2 ee2 ee3) 387 393 (* 9/68 (+ (* ee2 ee5) (* ee3 ee4)))))) 388 (+ (* 3 sigma)389 390 (expt an 3/2))))))394 (+ (* 3 sigma) 395 (/ (* power4 s) 396 (expt an 3/2)))))) 391 397 392 398 ;; Complete elliptic integral of the first kind. This can be computed … … 405 411 (t 406 412 (let ((precision (floatcontagion m))) 407 (carlsonrf (coerce 0 precision) ( 1 m) (coerce 1 precision))))))413 (carlsonrf 0 ( 1 m) 1))))) 408 414 409 415 ;; Elliptic integral of the first kind. This is computed using … … 418 424 Note for the complete elliptic integral, you can use elliptick" 419 425 (let* ((precision (floatcontagion x m)) 420 (x ( coercex precision))421 (m ( coercem precision)))426 (x (applycontagion x precision)) 427 (m (applycontagion m precision))) 422 428 (cond ((and (realp m) (realp x)) 423 429 (cond ((> m 1) … … 427 433 ;; 428 434 ;; with sin(theta) = sqrt(m)*sin(phi) 429 (/ (ellipticf ( cl:asin (* (sqrt m) (sin x))) (/ m))435 (/ (ellipticf (asin (* (sqrt m) (sin x))) (/ m)) 430 436 (sqrt m))) 431 437 ((< m 0) … … 447 453 ;; F(phi,1) = log(sec(phi)+tan(phi)) 448 454 ;; = log(tan(pi/4+pi/2)) 449 (log ( cl:tan (+ (/ x 2) (/ (floatpi x) 4)))))455 (log (tan (+ (/ x 2) (/ (floatpi x) 4))))) 450 456 ((minusp x) 451 457 ( (ellipticf ( x) m))) … … 464 470 (* ( 1 (* k sinx)) 465 471 (+ 1 (* k sinx))) 466 1 .0))))472 1)))) 467 473 ((< x (floatpi x)) 468 474 (+ (* 2 (elliptick m)) … … 487 493 E(phi, m) = integrate(sqrt(1m*sin(x)^2), x, 0, phi)" 488 494 (let* ((precision (floatcontagion phi m)) 489 (phi ( coercephi precision))490 (m ( coercem precision)))495 (phi (applycontagion phi precision)) 496 (m (applycontagion m precision))) 491 497 (cond ((= m 0) 492 498 ;; A&S 17.4.23 … … 502 508 (+ 1 (* k sinphi))))) 503 509 ( (* sinphi 504 (carlsonrf (* cosphi cosphi) y (coerce 1 precision)))510 (carlsonrf (* cosphi cosphi) y 1)) 505 511 (* (/ m 3) 506 512 (expt sinphi 3) 507 (carlsonrd (* cosphi cosphi) y (coerce 1 precision)))))))))513 (carlsonrd (* cosphi cosphi) y 1)))))))) 508 514 509 515 ;; Complete elliptic integral of second kind. … … 515 521 516 522 E(m) = integrate(sqrt(1m*sin(x)^2), x, 0, %pi/2)" 517 (let ((precision (floatcontagion m))) 518 (cond ((= m 0) 519 ;; A&S 17.4.23 520 (/ (floatpi m) 2)) 521 ((= m 1) 522 ;; A&S 17.4.25 523 (coerce 1 precision)) 524 (t 525 (let* ((k (sqrt m)) 526 (y (* ( 1 k) 527 (+ 1 k)))) 528 ( (carlsonrf 0.0 y 1.0) 529 (* (/ m 3) 530 (carlsonrd 0.0 y 1.0)))))))) 523 (cond ((= m 0) 524 ;; A&S 17.4.23 525 (/ (floatpi m) 2)) 526 ((= m 1) 527 ;; A&S 17.4.25 528 (float 1 m)) 529 (t 530 (let* ((k (sqrt m)) 531 (y (* ( 1 k) 532 (+ 1 k)))) 533 ( (carlsonrf 0 y 1) 534 (* (/ m 3) 535 (carlsonrd 0 y 1)))))))
Note: See TracChangeset
for help on using the changeset viewer.