archive/share/trash/ai.mc

   1 ;;;;;;;;;;;;;;;;;;; -*- Mode: Lisp; Package: Macsyma -*- ;;;;;;;;;;;;;;;;;;;
   2 ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
   3 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   4
   5 (macsyma-module ai)
   6
   7 (DEFGRAD $AI (X)
   8   (($DAI) X))
   9 (defgrad $dai (x)
  10   ((mtimes) x (($ai) x)))
  11 (defgrad $bi (x)
  12   (($dbi) x))
  13 (defgrad $dbi (x)
  14   ((mtimes) x (($bi) x)))
  15
  16 (declare (flonum (airy flonum))
  17          (flonum rz))
  18
  19 (defun airy (rz)
  20   (declare (flonum temp temp1 temp2))
  21   (cond ((< (abs rz) 3.0)
  22          (let ((temp (*$ rz rz rz)))
  23            (-$ (+$ .3550281
  24                    (*$ temp
  25                        (+$ .05917135
  26                            (*$ temp
  27                                (+$ 1.97237834e-3
  28                                    (*$ temp
  29                                        (+$ 2.73941433e-5
  30                                            (*$ temp
  31                                                (+$ 2.07531393e-7
  32                                                    (*$ temp
  33                                                        (+$ 9.8824472e-10
  34                                                            (*$ temp 3.22955787e-12))))))))))))
  35                (+$ (*$ rz
  36                        (+$ .2588194
  37                            (*$ temp
  38                                (+$ .0215682834
  39                                    (*$ temp
  40                                        (+$ 5.1353055e-4
  41                                            (*$ temp
  42                                                (+$ 5.70589507e-6
  43                                                    (*$ temp
  44                                                        (+$ 3.65762505e-8
  45                                                            (*$ temp
  46                                                                (+$ 1.52401045e-10
  47                                                                    (*$ temp 4.45617086e-13)))))))))))))))))
  48         ((> 0.0 rz)
  49          (let ((temp1 (*$ .66666667 (expt (*$ -1.0 rz) 1.5)))
  50                (temp2 (//$ -2.25 (*$ rz rz rz))))
  51            (*$ (//$ 1.0 (sqrt (*$ 3.1415926535 (sqrt (*$ -1.0 rz)))))
  52                (-$ (*$ (sin (+$ temp1 .78539816))
  53                        (+$ 1.0
  54                            (*$ temp2
  55                                (+$ -.037133488
  56                                    (*$ temp2 .0576491905)))))
  57                    (*$ (//$ (cos (+$ temp1 .78539816)) temp1)
  58                        (+$ .069444444
  59                            (*$ temp2
  60                                (+$ -.037993195
  61                                    (*$ temp2 .116099063)))))))))
  62         (t
  63          (let ((temp2 (//$ 1.5 (expt rz 1.5))))
  64            (*$ (//$ .5 (sqrt (*$ 3.1415926535 (sqrt rz))))
  65                (exp (*$ -1.0 (*$ .6666667 (expt rz 1.5))))
  66                (+$ 1.0
  67                    (*$ temp2
  68                        (+$ -.069444444
  69                            (*$ temp2
  70                                (+$ .037133488
  71                                    (*$ temp2
  72                                        (+$ -.037993195
  73                                            (*$ temp2
  74                                                (+$ .0576491905
  75                                                    (*$ temp2
  76                                                        (+$ -.116099063
  77                                                            (*$ temp2 .2915914)))))))))))))))))
  78
  79 (defun $ai ($arg)
  80   (cond ((numberp $arg) (airy (float $arg)))
  81         (t (list '($ai simp) $arg))))
  82
  83 (declare (flonum (bi flonum)))
  84
  85 (defun bairy (rz)
  86   (declare (flonum rz temp temp1 temp2))
  87   (cond ((< (abs rz) 3.0)
  88          (let ((temp (*$ rz rz rz)))
  89            (+$ (+$ .61492671
  90                    (*$ temp
  91                        (+$ .102487785
  92                            (*$ temp
  93                                (+$ 3.4162595e-3
  94                                    (*$ temp
  95                                        (+$ 4.74480486e-5
  96                                            (*$ temp
  97                                                (+$ 3.59454915e-7
  98                                                    (*$ temp
  99                                                        (+$ 1.71169007e-9
 100                                                            (*$ temp 5.59375834e-12))))))))))))
 101                (*$ rz
 102                    (+$ .44828835
 103                        (*$ temp
 104                            (+$ .0373573625
 105                                (*$ temp
 106                                    (+$ 8.8946102e-4
 107                                        (*$ temp
 108                                            (+$ 9.8829001e-6
 109                                                (*$ temp
 110                                                    (+$ 6.3351924e-8
 111                                                        (*$ temp
 112                                                            (+$ 2.6396635e-10
 113                                                                (*$ temp 7.71831435e-13))))))))))))))))
 114         ((> 0.0 rz)
 115          (let ((temp1 (*$ .66666667 (expt (*$ -1.0 rz) 1.5)))
 116                (temp2 (//$ -2.25 (*$ rz rz rz))))
 117            (*$ (//$ 1.0 (sqrt (*$ 3.1415926535 (sqrt (*$ -1.0 rz)))))
 118                (+$
 119                 (*$ (cos (+$ temp1 .78539816))
 120                     (+$ 1.0
 121                         (*$ temp2
 122                             (+$ -.037133488
 123                                 (*$ .0576491905 temp2)))))
 124                 (*$ (//$ (sin (+$ temp1 .78539816)) temp1)
 125                     (+$ .069444444
 126                         (*$ temp2
 127                             (+$ -.037993195
 128                                 (*$ .116099063 temp2)))))))))
 129         (t
 130          (let ((temp2 (//$ 1.5 (expt rz 1.5))))
 131            (*$ (//$ (exp (*$ .66666667 (expt rz 1.5)))
 132                     (sqrt (*$ 3.1415926535 (sqrt rz))))
 133                (+$ 1.0
 134                    (*$ temp2
 135                        (+$ .069444444
 136                            (*$ temp2
 137                                (+$ .037133488
 138                                    (*$ temp2
 139                                        (+$ .037993195
 140                                            (*$ temp2
 141                                                (+$ .0576491905
 142                                                    (*$ temp2
 143                                                        (+$ .116099063
 144                                                            (*$ temp2 .2915914)))))))))))))))))
 145
 146 (defun $bi ($arg)
 147   (cond ((numberp $arg) (bairy (float $arg)))
 148         (t (list '($bi simp) $arg))))
 149
 150
 151 (declare (flonum (dairy flonum)))
 152
 153 (defun dairy (rz)
 154   (declare (flonum rz temp temp1 temp2))
 155   (cond ((< (abs rz) 3.0)
 156          (let ((temp (*$ rz rz rz)))
 157            (-$ (*$ rz
 158                    rz
 159                    (+$ .17751405
 160                        (*$ temp
 161                            (+$ .01183427
 162                                (*$ temp
 163                                    (+$ 2.4654294e-4
 164                                        (*$ temp
 165                                            (+$ 2.49037668e-6
 166                                                (*$ temp
 167                                                    (+$ 1.48236707e-8
 168                                                        (*$ temp 5.8132042e-11)))))))))))
 169                (+$ .2588194
 170                    (*$ temp
 171                        (+$ .086273134
 172                            (*$ temp
 173                                (+$ 3.5947139e-3
 174                                    (*$ temp
 175                                        (+$ 5.70589507e-5
 176                                            (*$ temp
 177                                                (+$ 4.7549126e-7
 178                                                    (*$ temp
 179                                                        (+$ 2.43841672e-9
 180                                                            (*$ temp 8.4667248e-12)))))))))))))))
 181         ((> 0.0 rz)
 182          (let ((temp1 (*$ .66666667 (expt (*$ -1.0 rz) 1.5)))
 183                (temp2 (//$ -2.25 (*$ rz rz rz))))
 184            (*$ -1.0
 185                (sqrt (//$ (sqrt (*$ -1.0 rz)) 3.1415926535))
 186                (+$ (*$ (cos (+$ temp1 .78539816))
 187                        (+$ 1.0
 188                            (*$ temp2
 189                                (+$ .043885031
 190                                    (*$ temp2
 191                                        (+$ -.062662164
 192                                            (*$ temp2 3.08253244)))))))
 193                    (*$ (//$ (sin (+$ temp1 .78539816)) temp1)
 194                        (+$ -.097222222
 195                            (*$ temp2
 196                                (+$ .042462831
 197                                    (*$ temp2
 198                                        (+$ -.124105896
 199                                            (*$ temp2 .92047998)))))))))))
 200         (t
 201          (let ((temp2 (//$ 1.5 (expt rz 1.5))))
 202            (*$ -.5
 203                (sqrt (//$ (sqrt rz) 3.1415926535))
 204                (exp (*$ -.66666667 (expt rz 1.5)))
 205                (+$ 1.0
 206                    (*$ temp2
 207                        (+$ .097222222
 208                            (*$ temp2
 209                                (+$ -.043885031
 210                                    (*$ temp2
 211                                        (+$ .042462831
 212                                            (*$ temp2
 213                                                (+$ -.062662164
 214                                                    (*$ temp2
 215                                                        (+$ .124105896
 216                                                            (*$ temp2 -3.08253244)))))))))))))))))
 217
 218 (defun $dai ($arg)
 219   (cond ((numberp $arg) (dairy (float $arg)))
 220         (t (list '($dai simp) $arg))))
 221
 222 (declare (flonum (dbairy flonum)))
 223
 224 (defun dbairy (rz)
 225   (declare (flonum rz temp temp1 temp2))
 226   (cond ((< (abs rz) 3.0)
 227          (let ((temp (*$ rz rz rz)))
 228            (+$ (*$ rz
 229                    rz
 230                    (+$ .307463355
 231                        (*$ temp
 232                            (+$ .020497557
 233                                (*$ temp
 234                                    (+$ 4.2703244e-4
 235                                        (*$ temp
 236                                            (+$ 4.313459e-6
 237                                                (*$ temp
 238                                                    (+$ 2.5675351e-8
 239                                                        (*$ temp 1.00687651e-10)))))))))))
 240                (+$ .44828835
 241                    (*$ temp
 242                        (+$ .14942945
 243                            (*$ temp
 244                                (+$ 6.2262271e-3
 245                                    (*$ temp
 246                                        (+$ 9.8829001e-5
 247                                            (*$ temp
 248                                                (+$ 8.2357502e-7
 249                                                    (*$ temp
 250                                                        (+$ 4.22346157e-9
 251                                                            (*$ temp 1.46647972e-11)))))))))))))))
 252         ((> 0.0 rz)
 253          (let ((temp1 (*$ .66666667 (expt (*$ -1.0 rz) 1.5)))
 254                (temp2 (//$ -2.25 (*$ rz rz rz))))
 255            (*$ (sqrt (//$ (sqrt (*$ -1.0 rz)) 3.14159265353))
 256                (-$ (*$ (sin (+$ temp1 .78539816))
 257                        (+$ 1.0
 258                            (*$ temp2
 259                                (+$ .043885031
 260                                    (*$ temp2
 261                                        (+$ -.062662164
 262                                            (*$ temp2 3.08253244)))))))
 263                    (*$ (//$ (cos (+$ temp1 .78539816)) temp1)
 264                        (+$ -.097222222
 265                            (*$ temp2
 266                                (+$ .042462831
 267                                    (*$ temp2
 268                                        (+$ -.124105896
 269                                            (*$ temp2 .92037998)))))))))))
 270         (t
 271          (let ((temp2 (//$ 1.5 (expt rz 1.5))))
 272            (*$ (sqrt (//$ (sqrt rz) 3.1415926535))
 273                (exp (*$ .66666667 (expt rz 1.5)))
 274                (+$ 1.0
 275                    (*$ temp2
 276                        (+$ -.097222222
 277                            (*$ temp2
 278                                (+$ -.043885031
 279                                    (*$ temp2
 280                                        (+$ -.042462831
 281                                            (*$ temp2
 282                                                (+$ -.062662164
 283                                                    (*$ temp2
 284                                                        (+$ -.124105896
 285                                                            (*$ temp2 -3.08253244)))))))))))))))))
 286
 287 (defun $dbi ($arg)
 288   (cond ((numberp $arg) (dbairy (float $arg)))
 289         (t (list '($dbi simp) $arg))))