laplace.cpp

   1 // Laplace transform
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 eval_laplace(void)
   8 {
   9         push(cadr(p1));
  10         eval();
  11         push(symbol(SYMBOL_T));
  12         laplace();
  13 }
  14
  15 #define F p3
  16 #define T p4
  17 #define A p5
  18
  19 void
  20 laplace(void)
  21 {
  22         int h;
  23         save();
  24
  25         T = pop();
  26         F = pop();
  27
  28         // L[f + g] = L[f] + L[g]
  29
  30         if (car(F) == symbol(ADD)) {
  31                 p1 = cdr(F);
  32                 h = tos;
  33                 while (iscons(p1)) {
  34                         push(car(p1));
  35                         push(T);
  36                         laplace();
  37                         p1 = cdr(p1);
  38                 }
  39                 add_all(tos - h);
  40                 restore();
  41                 return;
  42         }
  43
  44         // L[Af] = A L[f]
  45
  46         if (car(F) == symbol(MULTIPLY)) {
  47                 push(F);
  48                 push(T);
  49                 partition();
  50                 F = pop();
  51                 A = pop();
  52                 laplace_main();
  53                 push(A);
  54                 multiply();
  55         } else
  56                 laplace_main();
  57
  58         restore();
  59 }
  60
  61 void
  62 laplace_main(void)
  63 {
  64         int n;
  65
  66         // L[t] = 1 / s^2
  67
  68         if (F == symbol(SYMBOL_T)) {
  69                 push_symbol(SYMBOL_S);
  70                 push_integer(-2);
  71                 power();
  72                 return;
  73         }
  74
  75         // L[t^n] = n! / s^(n+1)
  76
  77         if (car(F) == symbol(POWER) && cadr(F) == T) {
  78                 push(caddr(F));
  79                 n = pop_integer();
  80                 if (n > 0) {
  81                         push_integer(n);
  82                         factorial();
  83                         push_symbol(SYMBOL_S);
  84                         push_integer(n + 1);
  85                         power();
  86                         divide();
  87                         return;
  88                 }
  89         }
  90
  91         stop("laplace: cannot solve");
  92 }
  93
  94 #if SELFTEST
  95
  96 static char *s[] = {
  97
  98 // float ok?
  99
 100         "laplace(3t^2.0)",
 101         "6/(s^3)",
 102 };
 103
 104 void
 105 test_laplace(void)
 106 {
 107         test(__FILE__, s, sizeof s / sizeof (char *));
 108 }
 109
 110 #endif