src/plasma.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8
   9 (in-package :maxima)
  10
  11 ;; Plasma Dispersion Function, NZETA(Z).
  12 ;;
  13 ;; This function is related to the complex error function by
  14 ;;
  15 ;;    NZETA(Z) = %I*SQRT(%PI)*EXP(-Z^2)*(1-ERF(-%I*Z))
  16 ;;
  17 ;; NZETA(Z)  returns the complex value of the Plasma Dispersion Function
  18 ;; for complex Z.
  19 ;;
  20 ;; NZETAR(Z)  returns REALPART(NZETA(Z)).
  21 ;;
  22 ;; NZETAI(Z)  returns IMAGPART(NZETA(Z)).
  23
  24 (defun z-function (x y)
  25   (let ((xs (if (> 0.0 x) -1.0 1.0))
  26         (ys (if (> 0.0 y) -1.0 1.0))
  27         (capn 0) (nu 0)
  28         (bool nil)
  29         (h 0.0) (h2 0.0) (lamb 0.0) (r1 0.0) (r2 0.0) (s 0.0)
  30         (s1 0.0) (s2 0.0) (t1 0.0) (t2 0.0) (c 0.0)
  31         (re 0.0) (im 0.0))
  32     (setq x (abs x) y (abs y))
  33     (cond ((and (> 4.29 y) (> 5.33 x))
  34            (setq s (* (1+ (* -0.23310023 y))
  35                       (sqrt (1+ (* -0.035198873 x x)))))
  36            (setq h (* 1.6 s) h2 (* 2.0 h) capn (+ 6 (floor (* 23.0 s))))
  37            (setq nu (+ 9 (floor (* 21.0 s)))))
  38           (t (setq h 0.0) (setq capn 0) (setq nu 8)))
  39     (when (> h 0.0) (setq lamb (expt h2 capn)))
  40     (setq bool (or (zerop h) (zerop lamb)))
  41     (do ((n nu (1- n)))
  42         ((> 0 n))
  43       (setq t1 (+ h (* (float (1+ n)) r1) y))
  44       (setq t2 (- x (* (float (1+ n)) r2)))
  45       (setq c (/ 0.5 (+ (* t1 t1) (* t2 t2))))
  46       (setq r1 (* c t1) r2 (* c t2))
  47       (cond ((and (> h 0.0) (not (< capn n)))
  48              (setq t1 (+ s1 lamb) s1 (- (* r1 t1) (* r2 s2)))
  49              (setq s2 (+ (* r1 s2) (* r2 t1)) lamb (/ lamb h2)))))
  50     (setq im (if (zerop y)
  51                  (* 1.77245384 (exp (- (* x x))))
  52                  (* 2.0 (if bool r1 s1))))
  53     (setq re (* -2.0 (if bool r2 s2)))
  54     (cond ((> ys 0.0) (setq re (* re xs)))
  55           (t (setq r1 (* 3.5449077 (exp (- (* y y) (* x x)))))
  56              (setq r2 (* 2.0 x y))
  57              (setq re (* (- re (* r1 (sin r2))) xs))
  58              (setq im (- (* r1 (cos r2)) im))))
  59     `((mlist simp) ,re ,im)))
  60
  61 (defmfun $nzeta (z)
  62   (let ((x ($realpart z))
  63         (y ($imagpart z)))
  64      (if (and (numberp x) (numberp y))
  65          (let ((w (z-function (float x) (float y))))
  66            (simplify `((mplus) ,(second w) ,(simplify `((mtimes) $%i ,(third w))))))
  67          `(($nzeta simp) ,z))))
  68
  69 (defmfun $nzetar (z)
  70   (let ((x ($realpart z))
  71         (y ($imagpart z)))
  72      (if (and (numberp x) (numberp y))
  73          (second (z-function (float x) (float y)))
  74          `(($nzetar simp) ,z))))
  75
  76 (defmfun $nzetai (z)
  77   (let ((x ($realpart z))
  78         (y ($imagpart z)))
  79      (if (and (numberp x) (numberp y))
  80          (third (z-function (float x) (float y)))
  81          `(($nzetai simp) ,z))))