(*:Version: Mathematica 1.2 *)

(*:Name: Calculus`Laplace` *)

(*:Title: Laplace *)

(*:Author: Eran Yehudai *)

(*:Summery:

Implements Laplace and Inverse Laplace transforms.

*)

(*:Keywords:

laplace;
laplace transform;
differential equations;

*)

(*:Requirements: None *)

(*:Sources:
Fritz Oberhettinger and Larry Badii (1973), "Tables of Laplace
Transforms", New-York: Springer-Verlag.

*)

BeginPackage["Calculus`Laplace`"]
Laplace::usage =
  "Laplace[expr, t, s, opts] gives a function which is the Laplace
transform of expr. It is defined by Laplace[expr, t, s] =
Integrate[Exp[-s t] expr, {t, 0, Infinity}].
The option InfinityLimit determines how to treat the limit t->Infinity. "
Options[Laplace] = {InfinityLimit->Automatic}
SetAttributes[Laplace, Listable]
InverseLaplace::usage =
  "InverseLaplace[expr, s, t, opts] gives a function,
the Laplace transform of which is expr.
The option InfinityLimit determines how to treat the limit t->Infinity. "
Options[InverseLaplace] = {InfinityLimit->Automatic}
SetAttributes[InverseLaplace, Listable]
InfinityLimit::usage =
  "InfinityLimit is an option of Laplace. It determines how to treat
the limit t->Infinity of Integrate[Exp[-s t] fun, {t, 0, Infinity}].
InfinityLimit->All uses the Mathematica Limit function.
InfinityLimit->Hold returns the limit in the form
Limit[expr, Hold[t->Infinity].
InfinityLimit->Automatic (default) is the same as InfinityLimit->All
unless either error messages are produced or the final expressions
contains Infinity. In these cases, the result is the same as with
InfinityLimit->Hold. InfinityLimit->func applies func to the expression.
InfinityLimit->Limit[#1, #2->Infinity]& is the same as InfinityLimit->All.
In particular, InfinityLimit->(0&) often gives the right result.
A notable exception is when the original function contains terms
of the form Exp[-t^2]."
Delta::usage =
  "Delta[expr] is Chroniker's delta, defined by
Integrate[Delta[t] f[t], {t, t1, t2}] = If[t1 <= 0 <= t2, f[0], 0]"
DeltaRules::usage =
  "DeltaRules are the rules that implement a primitive
form of integration over the Delta function. Use by
Integrate[__] //. DeltaRules"
ComposedFunctionQ::usage =
  "ComposedFunctionQ[expr, x] return a list of the smallest subexpressions
containing x in expr."
ThomsonBei::usage =
  "ThomsonBei[nu, z] is the Thomson bei function."
ThomsonBer::usage =
  "ThomsonBer[nu, z] is the Thomson ber function."
ParabolicCylinderD::usage =
  "ParabolicCylinderD[nu, z] is the D parabolic cylinder function."
DefiniteIntegrateZeroToInfinity::usage =
  "DefiniteIntegrateZeroToInfinity[fun, t, defopt] is the definite
integral from zero to infinity of f as a function of t, with the
limit to infinity treated according to defopt."
Options[DefiniteIntegrateZeroToInfinity] := {InfinityLimit->Automatic}
(************************************************************************)
Begin["`Private`"]
(************************************************************************)
(* ==== Thomson Functions ============================================= *)
ThomsonBei[n_, z_] :=
  (BesselJ[n, z Exp[3 Pi I/4]] - BesselJ[n, z Exp[-3 Pi I/4]]) / (2I)
ThomsonBei[z_] := ThomsonBei[0, z]
ThomsonBer[n_, z_] :=
  (BesselJ[n, z Exp[3 Pi I/4]] + BesselJ[n, z Exp[-3 Pi I/4]]) / 2
ThomsonBer[z_] := ThomsonBer[0, z]
(* ==== ParabolicCylinderD ============================================ *)
ParabolicCylinderD[-1,z_] := Sqrt[Pi/2]Exp[z^2/4](1-Erf[z/Sqrt[2]])
ParabolicCylinderD[-1/2,z_] := Sqrt[z/(2Pi)] BesselK[z^2/4]
ParabolicCylinderD[p_, z_] := 2^(p/2)Exp[-z^2/4] (
  Sqrt[Pi]Hypergeometric1F1[-p/2,1/2,z^2/2]/Gamma[(1-p)/2] -
  Sqrt[2Pi] z Hypergeometric1F1[(1-p)/2,3/2,z^2/2]/Gamma[-p/2])
(* ==== Delta Function ================================================ *)
DeltaRules =
  {
    Integrate[c_. + expr_. Delta[a_. t_ + b_.], t_] :>
      Integrate[c, t] + If[t > -b/a, expr/a, 0] /; FreeQ[{a, b}, t],
    Integrate[c_. + expr_. Delta[a_. t_ + b_.], {t_, t1_, t2_}] :>
      Integrate[c, {t, t1, t2}] +
        Simplify[If[t1 <= -b/a <= t2, Limit[expr, t-> -b/a]/a, 0]] /;
        FreeQ[{a, b}, t]
  }
Delta[x_] := 0 /; x != 0
(* ==== DefiniteIntegrateZeroToInfinity =============================== *)
DefiniteIntegrateZeroToInfinity[f1_ + f2_, t_, opt___] :=
  DefiniteIntegrateZeroToInfinity[f1, t, opt] +
    DefiniteIntegrateZeroToInfinity[f2, t, opt]
DefiniteIntegrateZeroToInfinity[expr_. Delta[a_. t_ + b_.], t_, ___] :=
  If[Positive[a b], 0, Limit[expr, t -> -b/a]]
DefiniteIntegrateZeroToInfinity[f_, t_, opt___] :=
  Block[
    {
      il = InfinityLimit /. {opt} /.
        Options[DefiniteIntegrateZeroToInfinity],
      g = Integrate[Collect[f, Delta], t] //. DeltaRules,
      g0,
      gi
    },
    If[
      FreeQ[g, Integrate],
      (
        Off[Limit::nlm, Limit::lim];
        g0 = Limit[g, t->0];
        holdForm = Limit[g, Hold[t->Infinity]];
        gi =
          Switch[
            il,
            Hold,
            holdForm,
            All,
            Limit[g, t->Infinity],
            Automatic,
            Block[
              {
                ggi =
                  (
                    Off[
                      Limit::nlm,
                      Limit::lim,
                      SeriesData::csa,
                      Infinity::indt
                    ];
                    Check[
                      Limit[
                        g //. {x_Plus t :> Expand[x t]},
                        t->Infinity
                      ] //.
                          {
                            Erf[_. -Infinity] -> -1,
                            Erf[_. Infinity]->1,
                            __ Exp[_. + -Infinity __] -> 0
                          },
                      holdForm
                    ]
                  )
              },
              If[
                Head[ggi] =!= Limit && !FreeQ[ggi, DirectedInfinity],
                holdForm,
                ggi
              ]
            ],
            _,
            il[g, t]
          ];
        gi-g0
      ),
      Integrate[Collect[f, Delta], {t, 0, Infinity}] //. DeltaRules
    ]
  ]
(* ==== ComposedFunctionQ ============================================ *)
ComposedFunctionQ[f_, s_Symbol] :=
  If[
    FreeQ[f, s] || f===s,
    False,
    (Drop[#,-1]& /@ Position[f, s]) /. {aaaaa___Integer} :> f[[aaaaa]]
  ]
(************************************************************************)
(*                        Laplace Transform                             *)(************************************************************************)

(* ==== Numeric Value ================================================= *)
N[Laplace[f_, t_, s_?NumberQ, ___]] := NIntegrate[Exp[-s t]f,{t,0,Infinity}]
N[Laplace[f_, t_, s_?NumberQ, ___], d_] :=
  NIntegrate[Exp[-s t]f,{t,0,Infinity}, WorkingPrecision->d]
(* ==== Using Apart =================================================== *)
Laplace[f:Power[x_, _?Negative], t_, s_, opt___] :=
  Block[{ff = Apart[f, t]}, Laplace[ff, t, s, opt] /; ff =!= f] /; !FreeQ[x, t]
Laplace[f:Times[__, Power[x_, _?Negative]], t_, s_, opt___] :=
  Block[{ff = Apart[f, t]}, Laplace[ff, t, s, opt] /; ff =!= f] /; !FreeQ[x, t]
(* ==== Linearity ===================================================== *)
Laplace[f_ + g_, t_, s_, opt___] :=
  Laplace[f, t, s, opt] + Laplace[g, t, s, opt]
Laplace[c_ f_, t_, s_, opt___] := c Laplace[f, t, s, opt] /; FreeQ[c, t]
Laplace[x_Plus f_, t_, s_, opt___] := Laplace[Expand[x f], t, s, opt]
(* ==== Products Involving Powers ===================================== *)
Laplace[a_ t_, t_, s_, opt___] := -D[Laplace[a, t, s, opt], s]
Laplace[a_ t_^(n_Integer), t_, s_, opt___] :=
  (-1)^n D[Laplace[a, t, s, opt], {s, n}] /; (FreeQ[n, t] && n > 0)
Laplace[t_^n_., t_, s_, ___] := Gamma[n+1] s^(-n-1) /; FreeQ[n, t]
Laplace[a_/t_, t_, s_, opt___] :=
  Block[
    {v = Unique["t"]},
    Integrate[Laplace[a, t, v, opt],{v, s, Infinity}]
  ]
(* ==== Other Rational Functions ====================================== *)
Laplace[1/(a_. t_ + b_), t_, s_, ___] :=
  -Exp[b s/a]ExpIntegralEi[-b s/a]/a /; FreeQ[{a,b},t]
Laplace[1 /(d_. t_^2+a_?Negative b_.), t_, s_, ___] :=
  Block[
    {c = Sqrt[-a b/d]},
    ExpIntegralEi[-c s]Exp[c s]/(2c d) - ExpIntegralEi[c s]Exp[-c s]/(2c d) /;
    FreeQ[{d, c}, t]
  ]
Laplace[(a_. t_+b_.)/(f_. t_^2+c_?Negative d_.), t_, s_, ___] :=
  Block[
    {e = Sqrt[-c d/f]},
    -ExpIntegralEi[-e s]Exp[e s](a-b/e)/(2f) -
      ExpIntegralEi[e s]Exp[-e s](a+b/e)/(2f) /;
    FreeQ[{a, b, e, f}, t]
  ]
Laplace[1/(a_. t_^2+c_), t_, s_, ___] :=
  CosIntegral[Sqrt[c/a] s] Sin[Sqrt[c/a]s]/(Sqrt[a c]) +
    (Pi/2-SinIntegral[Sqrt[c/a] s]) Cos[Sqrt[c/a]s]/(Sqrt[a c]) /;
  FreeQ[{a, c}, t]
Laplace[(a_. t_+b_.)/(d_. t_^2+c_), t_, s_, ___] :=
  -CosIntegral[Sqrt[c/d] s]((a/d) Cos[Sqrt[c/d]s] -
    (b/Sqrt[c d])Sin[Sqrt[c/d]s]) +
    (Pi/2-SinIntegral[Sqrt[c/d] s]) *
    ((a/d) Sin[Sqrt[c/d]s]+(b/Sqrt[c d])Cos[Sqrt[c/d]s]) /;
  FreeQ[{a, b, c, d}, t]
(* ==== Exponentials================================================== *)
Laplace[a_. Power[E, (b_.+ c_. t_)], t_, s_, opt___] :=
  Laplace[a Exp[b], t, s-c, opt] /; FreeQ[{b, c}, t]
Laplace[Exp[a_?Negative b_. / t_], t_, s_, ___] :=
  Block[
    {c= -a b},
    2Sqrt[c/s] BesselK[1, 2Sqrt[c s]] /;
    FreeQ[c, t]
  ]
Laplace[t_^n_. Exp[a_?Negative b_. / t_], t_, s_, ___] :=
  Block[
    {nu=n+1, c= -a b},
    2(c/s)^((n+1)/2) BesselK[n+1, 2 Sqrt[c s]] /;
    FreeQ[{n, c}, t]
  ]
Laplace[Exp[a_?Negative b_. Sqrt[t_]], t_, s_, ___] :=
  Block[
    {c = -a b},
    s^(-1)Exp[c^2/(8s)]ParabolicCylinderD[-2, c/Sqrt[2s]] /;
    FreeQ[c, t]
  ]
Laplace[t_^n_. Exp[a_?Negative b_. Sqrt[t_]], t_, s_, ___] :=
  Block[
    {c = -a b},
    2(2s)^(-n-1)Gamma[2(n+1)] Exp[c^2/(8s)] *
      ParabolicCylinderD[-2(n+1), c/Sqrt[2s]] /;
    FreeQ[{n, c}, t]
  ]
Laplace[Exp[a_?Negative b_. t_^2], t_, s_, ___] :=
  Block[
    {c = -a b},
    (2 c)^(-1/2)Exp[s^2/(8c)] *
      ParabolicCylinderD[-1, s/Sqrt[2c]] /;
    FreeQ[c, t]
  ]
Laplace[t_^n_. Exp[a_?Negative b_. t_^2], t_, s_, ___] :=
  Block[
    {c = -a b},
    (2 c)^(-(n+1)/2)Gamma[n+1] Exp[s^2/(8c)] *
      ParabolicCylinderD[-(n+1), s/Sqrt[2c]] /;
    FreeQ[{n, c}, t]
  ]
Laplace[Exp[a_?Negative b_. Exp[t_]], t_, s_, ___] :=
  Block[{c = -a b}, c^s Gamma[-s, c] /; FreeQ[c,t]]
(* ==== Logarithms ==================================================== *)
Laplace[Log[t_], t_, s_, ___] := -(EulerGamma+Log[s])/s
Laplace[Log[a_+b_. t_], t_, s_, ___] :=
  (Log[a]-Exp[a s/b]ExpIntegralEi[-a s/b])/s /; FreeQ[{a,b},t]
Laplace[Log[t_]^2, t_, s_, ___] := (Pi^2/6+(EulerGamma+Log[s])^2)/s
Laplace[Log[a_+t_^2], t_, s_, ___] :=
  Block[
    {a = Sqrt[aa]},
    (2/s)(Log[a]-CosIntegral[a s]Cos[a X]-
      (SinIntegral[a s]-Pi/2)Sin[a s]) /;
    FreeQ[a,t]
  ]
Laplace[t_^n_. Log[t_], t_, s_, ___] :=
  Gamma[n+1]s^(-n-1)(PolyGamma[n+1]-Log[s]) /; FreeQ[n,t]
Laplace[t_^n_. Log[t_]^2, t_, s_, ___] :=
  Gamma[n+1]s^(-n-1)((PolyGamma[n+1]-Log[s])^2+PolyGamma[1,n+1]) /;
    FreeQ[n,t]
(* ==== Trigonometric Functions ======================================= *)
Laplace[t_^n_. Sin[a_. t_], t_, s_, ___] :=
  Gamma[n+1](s^2+a^2)^(-(n+1)/2) Sin[(n+1)ArcTan[a/s]] /; FreeQ[{n,a},t]
Laplace[Sin[a_. t_], t_, s_, ___] := a/(a^2+s^2)
Laplace[t_^n_. Cos[a_. t_], t_, s_, ___] :=
  Gamma[n+1](s^2+a^2)^(-(n+1)/2) Cos[(n+1)ArcTan[a/s]] /; FreeQ[{n,a},t]
Laplace[Cos[a_. t_], t_, s_, ___] := s/(a^2+s^2)
(* ==== Inverse Trigonometric Functions =============================== *)
Laplace[ArcTan[a_. t_], t_, s_, ___] :=
  (CosIntegral[s/a]Sin[s/a]- (SinIntegral[s/a]-Pi/2)Cos[s/a])/s /; FreeQ[a,t]
Laplace[ArcCot[a_. t_], t_, s_, ___] :=
  ((Pi/2-CosIntegral[s/a])Sin[s/a] - (SinIntegral[s/a]-Pi/2)Cos[s/a])/s /;
    FreeQ[a,t]
(* ==== Hyperbolic Functions ========================================== *)
Laplace[Tanh[a_. t_], t_, s_, ___] :=
  (PolyGamma[1/2+s/(4a)]-PolyGamma[s/(4a)])/(2a)- 1/s /; FreeQ[a,t]
(* ==== Derivatives and Integrals ===================================== *)
Laplace[Derivative[1][f_][t_], t_, s_, opt___] :=
  s Laplace[f[t], t, s, opt] - f[0]
Laplace[Derivative[n_Integer][f_][t_], t_, s_, opt___] :=
  s^n Laplace[f[t], t, s, opt] - (Sum[s^i D[f[t], {t, n-1-i}] ,
    			{i, 0, n-1}] /. t->0)Laplace[Integrate[f_, t_], t_, s_, opt___] :=
  Laplace[f, t, s, opt] / s
(* ==== InverseLaplace ================================================ *)
Laplace[InverseLaplace[f_, s_, t_, ___], t_, s1_, ___] :=
  f /. s->s1
(* ==== Delta ========================================================= *)
Laplace[expr_. Delta[a_. t_ + b_.], t_, s_, ___] :=
  If[a b > 0, 0, Limit[expr Exp[-s t], t -> -b/a]]
(* ==== Using the Definition ========================================== *)
Laplace[f_, t_Symbol, s_, opt___] :=
  Block[
    {
      g =
        DefiniteIntegrateZeroToInfinity[
          Exp[-s t]f,
          t,
          InfinityLimit->(InfinityLimit /. {opt} /. Options[Laplace])
        ]
    },
    g /;
    FreeQ[g, Integrate]
  ]
(* ==== Composed Function Case ======================================== *)
Laplace[f_, t_, s_, opt___] :=
  Block[
    {g = ComposedFunctionQ[f, t]},
    Switch[
      g[[1]],
      _Times,
      Block[{a = g[[1]]/t}, Laplace[f /. g[[1]]->t, t, s/a, opt] / a]
    ] /;
    MatchQ[g[[1]], _Times] && SameQ @@ g
  ]


(************************************************************************)
(*                      Inverse Laplace Transform                       *)(************************************************************************)

InverseLaplace[f:Power[x_, _?Negative], s_, t_, opt___] :=
  Block[{ff = Apart[f, s]}, InverseLaplace[ff, s, t, opt] /; ff =!= f] /;
    !FreeQ[x, s]
InverseLaplace[f:Times[__, Power[x_, _?Negative]], s_, t_, opt___] :=
  Block[{ff = Apart[f, s]}, InverseLaplace[ff, s, t, opt] /; ff =!= f] /;
    !FreeQ[x, s]
(* ==== Linearity ===================================================== *)
InverseLaplace[c_, s_, t_, opt___] := c Delta[t] /; FreeQ[c, s]
InverseLaplace[c_ f_, s_, t_, opt___] := c InverseLaplace[f, s, t, opt] /;
  FreeQ[c, s]
InverseLaplace[x_Plus, s_, t_, opt___] := InverseLaplace[#, s, t, opt]& /@ x
InverseLaplace[a_Plus b_, s_, t_, opt___] :=
  InverseLaplace[#, s, t, opt]& /@ Expand[a b]
InverseLaplace[s_^(n1_?Negative n2_.), s_, t_, opt___] :=
  1/(t^(n1 n2+1) Gamma[-n1 n2]) /; FreeQ[{n1, n2}, s]
(* ==== DI ============================================================ *)
DI[f_, _, 0] := f
DI[f_, x_, n_Integer] :=
  Block[{v = Unique["t"]},
    Integrate[
      Collect[DI[f /. x->v, v, n-1], Delta],
      {v, 0, x}
    ] //. DeltaRules
  ] /; n > 0
(* ==== Other Rational Functions ====================================== *)
InverseLaplace[1/s_, s_, __] := 1
InverseLaplace[1/(s_+a_), s_, t_, opt___] := Exp[-a t] /; FreeQ[a, s]
InverseLaplace[(s_+a_.)^(n_Integer), s_, t_, opt___] :=
  t^(-1-n)Exp[-a t]/(-n-1)! /; (FreeQ[a, s] && Negative[n])
SqrtRule = {
             Sqrt[x_?Negative] :> I Sqrt[-x],
             1/Sqrt[x_?Negative] :> -I / Sqrt[-x]
           }
ExpRules = {
             Exp[x_. + Complex[r_, i_] y_.] :>
               Expand[Exp[x + r y](Cos[i y] + I Sin[i y])],
             Sin[x_?Negative y_.] :> -Sin[-x y],
             Cos[x_?Negative y_.] :> Cos[-x y]
           }
InverseLaplace[1/(a_. s_^2 + b_. s_ + c_.), s_, t_, opt___] :=
  Block[
    {
      x,
      v = ((x /. Solve[a x^2+b x+c==0, x]) /. SqrtRule)
    },
    ExpandAll[
      If[
        SameQ @@ v,
        t Exp[v[[1]] t],
        a (Exp[v[[1]] t]-Exp[v[[2]] t])/(Sqrt[b^2-4 a c] /. SqrtRule)
      ] /. Exp[xx:Times[_Plus, __]] :> Exp[Expand[xx]] //. ExpRules
    ]
  ]
InverseLaplace[1/(s_^2+c_?Negative b_.), s_, t_, opt___] :=
  Block[{a = Sqrt[-b c]}, Sinh[a t]/a /; FreeQ[a, s]]
InverseLaplace[1/(s_^2+b_), s_, t_, opt___] :=
  Block[{a = Sqrt[b]}, Sin[a t]/a /; FreeQ[a, s]]
InverseLaplace[1/(a_. s_^n_Integer + b_), s_, t_, opt___] :=
  Block[
    {c = Power[b/a, 1/n]},
    (c^(1-n)/(a n))
      Sum[Block[{z = Exp[I Pi(2k-1)/n]}, z Exp[-c z t] //. ExpRules], {k, n}]
  ] /; n > 0
InverseLaplace[Sqrt[s_+b_.]/(s_+a_.), s_, t_, opt___] :=
  Exp[-b t]/Sqrt[Pi t] + Sqrt[b-a] Exp[-a t] Erf[Sqrt[(b-a)t]] /;
    FreeQ[{a, b}, s]
InverseLaplace[1/Sqrt[s_+a_.]/Sqrt[s_+b_.], s_, t_, opt___] :=
  Exp[-(a+b)t/2] BesselI[0, (a-b)t/2]
InverseLaplace[(s_+a_.)^(-1/2)(s_+b_.)^(-3/2), s_, t_, opt___] :=
  t Exp[-(a+b)t/2] (BesselI[0, (a-b)t/2] + BesselI[1, (a-b)t/2])
InverseLaplace[Sqrt[s_+a_.]/(s_+b_.)^(3/2), s_, t_, opt___] :=
  Exp[-(a+b)t/2]((1+(a-b)t)BesselI[0, (a-b)t/2] + (a-b)t BesselI[1, (a-b)t/2])
InverseLaplace[(a_. s_^2+b_. s_+c_.)^Rational[k_?Negative, 2],s_,t_,opt___] :=
    InverseLaplace[
      ((s + b / (2 a))^2 + (4 a c - b^2) / (4 a))^(k/2),
      s,
      t,
      opt
    ] / a
InverseLaplace[
           (d_. s_^2+b_?Negative c_.)^Rational[k_?Negative,2],s_,t_,opt___] :=
  Block[
    {n = -k/2-1/2, a = Sqrt[-b c/d]},
    t^n BesselI[n, a t]/(If[n==0, 1, (2n-1)!!] a^n d)
  ] /; FreeQ[{b, c, d, k}, s]
InverseLaplace[(d_. s_^2 + c_)^Rational[k_?Negative, 2], s_, t_, opt___] :=
  Block[
    {n = -k/2-1/2, a = Sqrt[c/d]},
    t^n BesselJ[n, a t]/(If[n==0, 1, (2n-1)!!] a^n d)
  ] /; FreeQ[{c, d, k}, s]
InverseLaplace[(a_. s_ + b_. Sqrt[c_. s_^2 + d_])^n_., s_, t_, opt___] :=
  (a^n) InverseLaplace[(s + (b Sqrt[c] / a) Sqrt[s^2 + d/c])^n, s, t, opt] /;
    FreeQ[{a, b, c, d, n}, s] && (a =!= 1 || c =!= 1)
InverseLaplace[(s_ - Sqrt[s_^2 + c1_?Negative c2_.])^n_., s_, t_, opt___] :=
  (n/t) Sqrt[-c1 c2]^n BesselI[n, Sqrt[-c1 c2] t] /; FreeQ[{c1, c2, n}, s]
InverseLaplace[(s_ - Sqrt[s_^2 + c_])^n_., s_, t_, opt___] :=
  (-1)^n (n/t) Sqrt[c]^n BesselJ[n, Sqrt[c] t] /; FreeQ[{c, n}, s]
InverseLaplace[(c_. s_^2+b1_?Negative b2_.)^(n1_?Negative n2_.),s_,t_,opt___] :=
  Block[
    {a = Sqrt[-b1 b2/c], n = -n1 n2},
    Sqrt[Pi] (2a)^(1/2-n) / Gamma[n] t^(n-1/2) BesselI[n-1/2, a t] / c^n
  ] /; FreeQ[{c, b1, b2, n1, n2}, s]
InverseLaplace[(c_. s_^2 + b_)^(n1_?Negative n2_.), s_, t_, opt___] :=
  Block[
    {a = Sqrt[b/c], n = -n1 n2},
    Sqrt[Pi] (2a)^(1/2-n) / Gamma[n] t^(n-1/2) BesselJ[n-1/2, a t] / c^n
  ] /; FreeQ[{c, b, n1, n2}, s]
(* ==== Exponantials ================================================== *)
InverseLaplace[Exp[a_. s_ + b_.], s_, t_, opt___] :=
  Exp[b] Delta[t+a] /; FreeQ[{a,b}, s]
InverseLaplace[Exp[a1_?Negative a2_. / s_ + b_.], s_, t_, opt___] :=
  Block[
    {a = -a1 a2},
    Expand[Exp[b] (Delta[t] - Sqrt[a/t] BesselJ[1, 2 Sqrt[a t]])]
  ] /; FreeQ[{a1, a2, b}, s]
InverseLaplace[Exp[a_. / s_ + b_.], s_, t_, opt___] :=
  Exp[b] (Sqrt[a/t] BesselI[1, 2 Sqrt[a t]] + Delta[t]) /; FreeQ[{a, b}, s]
InverseLaplace[
    s_^(n1_?Negative n2_.) Exp[a1_?Negative a2_./s_ + b_.], s_, t_, opt___] :=
  Block[
    {n = -n1 n2 - 1, a = -a1 a2},
    Exp[b] (t/a)^(n/2) BesselJ[n, 2Sqrt[a t]]
  ] /; FreeQ[{n1, n2, a1, a2, b}, s]
InverseLaplace[s_^(n1_?Negative n2_.) Exp[a_./s_ + b_.], s_, t_, opt___] :=
  Block[
    {n = -n1 n2 - 1},
    Exp[b] (t/a)^(n/2) BesselI[n, 2Sqrt[a t]]
  ] /; FreeQ[{n1, n2, a, b}, s]
InverseLaplace[Exp[a1_?Negative a2_. Sqrt[s_] + b_.], s_, t_, opt___] :=
  Block[{a = -a1 a2}, Exp[b - a^2/(4t)] a/(2 Sqrt[Pi t^3])] /;
    FreeQ[{a1, a2, b}, s]
InverseLaplace[s_^n_. Exp[a1_?Negative a2_. Sqrt[s_] + b_.], s_, t_, opt___] :=
  Block[
    {a = -a1 a2},
    Exp[b] 2^(-n-1/2) Pi^(-1/2) t^(-n-1) Exp[-a^2/(8t)] *
      ParabolicCylinderD[2n+1, a/Sqrt[2t]]
  ] /; FreeQ[{a1, a2, b, n}, s]
(* ==== Logarithmic Functions ========================================= *)
InverseLaplace[s_^(n1_?Negative n2_.) Log[s_], s_, t_, opt___] :=
  Block[{n = -n1 n2}, t^(n-1)(PolyGamma[n]-Log[t])/Gamma[n]] /;
    FreeQ[{n1, n2}, s]
InverseLaplace[Log[a_. s_ + b_] s_^(n_Integer), s_, t_, opt___] :=
  DI[Log[b] - ExpIntegralEi[-b t/a], t, -1-n] /; FreeQ[{a, b}, s] && n < 0
InverseLaplace[Log[e_. + (a_. s_ + b_.) / (c_. s_ + d_.)], s_, t_, opt___] :=
  Log[e + a/c] Delta[t] + (-Exp[(d e + b) t/(c e + a)] + Exp[d t/c])/t /;
    FreeQ[{a, b, c, d, e}, s]
InverseLaplace[Log[e_. + b_. / (c_. s_ + d_.)], s_, t_, opt___] :=
  Log[e] Delta[t] + (-Exp[(d e + b) t/(c e)] + Exp[d t/c])/t /;
    FreeQ[{b, c, d, e}, s]
InverseLaplace[s_^(n1_?Negative n2_.) Log[s_]^2, s_, t_, opt___] :=
  Block[
    {n = -n1 n2},
    t^(n-1)((PolyGamma[n]-Log[t])^2 - PolyGamma[1, n]) / Gamma[n]
  ] /; FreeQ[{n1, n2}, s]
ilLog[a_, b_, c_, d_, t_] :=
  (2/t)(Cos[Sqrt[-c^2/4+d]t]Exp[-c t/2]-Cos[Sqrt[-a^2/4+b]t]Exp[-a t/2]) //.
    {
      Cos[Complex[aa_, bb_] cc_.] :>
        Cos[aa cc]Cosh[bb cc] - I Sin[aa cc]Sinh[bb cc],
      Cos[Sqrt[a1_?Negative] a2_.] :> Cosh[Sqrt[Expand[-a1]] a2]
    }
InverseLaplace[Log[d_. (a_. s_^2+c_.)/s_^2], s_, t_, opt___] :=
  Log[a d]Delta[t] + ilLog[0, c/a, 0, 0, t] /; FreeQ[{a,c,d},s]
InverseLaplace[Log[(a_. s_^2+c_.)/(d_. s_^2+f_.)], s_, t_, opt___] :=
  Log[a/d]Delta[t] + ilLog[0, c/a, 0, f/d, t] /; FreeQ[{a,c,d,f},s]
InverseLaplace[Log[(a_. s_^2+c_.)/(d_. s_^2+e_. s_+f_.)], s_, t_, opt___] :=
  Log[a/d]Delta[t] + ilLog[0, c/a, e/d, f/d, t] /; FreeQ[{a,c,d,e,f},s]
InverseLaplace[Log[(a_. s_^2+b_. s_+c_.)/(d_. s_^2+f_.)], s_, t_, opt___] :=
  Log[a/d]Delta[t] + ilLog[b/a, c/a, 0, f/d, t] /; FreeQ[{a,b,c,d,f},s]
InverseLaplace[Log[(a_. s_^2+b_. s_+c_.)/(d_. s_^2+e_. s_+f_.)],s_,t_,opt___] :=
  Log[a/d]Delta[t] + ilLog[b/a, c/a, e/d, f/d, t] /; FreeQ[{a,b,c,d,e,f},s]
InverseLaplace[
          Log[expr:Plus[a_.,Times[b_.,Power[_,_?Negative]]]],s_,t_,opt___] :=
  Block[
    {
      expr1 =
        MapAt[
          Collect[#, s]&,
          ExpandNumerator[ExpandDenominator[Together[expr]]],
          {{2, 1}, {1}}
        ]
    },
    InverseLaplace[Log[expr1], s, t, opt] /;
    expr1 =!= expr
  ]
(* ==== Trigonometric Functions ======================================= *)
InverseLaplace[s_^(n1_?Negative n2_.) Sin[a_. / s_], s_, t_, opt___] :=
  Block[
    {n = -n1 n2},
    (t/a)^((n-1)/2) (Sin[3Pi n/4 + Pi/4] ThomsonBer[n-1, 2Sqrt[a t]] -
      Cos[3Pi n/4+Pi/4] ThomsonBei[n-1, 2Sqrt[a t]])
  ] /; FreeQ[{n1, n2, a}, s]
InverseLaplace[s_^(n1_?Negative n2_.) Cos[a_. / s_], s_, t_, opt___] :=
  Block[
    {n = -n1 n2},
    -(t/a)^((n-1)/2) (Cos[3Pi n/4 + Pi/4] ThomsonBer[n-1, 2Sqrt[a t]] -
      Sin[3Pi n/4+Pi/4] ThomsonBei[n-1, 2Sqrt[a t]])
  ] /; FreeQ[{n1, n2, a}, s]
InverseLaplace[ArcTan[a_./s_], s_, t_, opt___] := Sin[a t] / t /; FreeQ[a, s]
(* ==== Derivatives and Integrals ===================================== *)
InverseLaplace[Derivative[n_Integer][f_][s_], s_, t_, opt___] :=
  Block[
    {g = InverseLaplace[f[s], s, t, opt]},
    (-1)^n t^n g /;
    Head[g] =!= InverseLaplace
  ]
InverseLaplace[a_ s_, s_, t_, opt___] :=
  Block[
    {f = InverseLaplace[a, s, t, opt]},
    D[f, t] + Limit[f, t->0] /;
    Head[f] =!= InverseLaplace
  ]
InverseLaplace[a_ s_^n_, s_, t_, opt___] :=
  Block[
    {f = InverseLaplace[a, s, t, opt]},
    If[
      n > 0,
      D[f, {t, n}] + Sum[Limit[D[f, {t, k}], t->0] t^(n-k-1), {k, 0, n-1}],
      DI[f, t, -n]
    ] /; FreeQ[f, InverseLaplace]
  ] /; Head[n] === Integer
(* ==== Laplace ========================================================== *)
InverseLaplace[Laplace[f_, t_, s_, ___], s_, t1_, ___] :=
  f /. t->t1
(* ==== Composed Functions =============================================== *)
InverseLaplace[g1_ g2_, s_, t_, opt___] :=
  Block[
    {
      f1 = InverseLaplace[g1, s, t, opt],
      f2 = InverseLaplace[g2, s, t, opt],
      v = Unique["t"]
    },
    (Integrate[
      Collect[(f1 /. t->v)(f2 /. t->t-v), Delta],
      {v, 0, t}
    ]  //.  DeltaRules) /;
    FreeQ[{f1, f2}, InverseLaplace]
  ]
InverseLaplace[f_, s_, t_, opt___] :=
  Block[
    {
      gg =
        Block[
          {g = ComposedFunctionQ[f, s]},
          Switch[
            g[[1]],
            _Times,
            Block[
              {a = g[[1]]/s},
              InverseLaplace[f /. g[[1]]->s, s, t/a, opt] / a
            ],
            _Plus,
            Block[
              {a = g[[1]]-s},
              InverseLaplace[f /. g[[1]]->s, s, t, opt] Exp[-a t]
            ],
            _,
            Block[
              {ff = InverseLaplace[f /. s->s^2, s, t, opt], v = Unique["t"]},
              Expand[
                Block[
                  {
                    di =
                      DefiniteIntegrateZeroToInfinity[
                        v Exp[-v^2/(4 t)] (ff /. t->v),
                        v,
                        InfinityLimit->
                          (InfinityLimit /.
                            {opt} /.Options[InverseLaplace])
                      ]
                  },
                  If[
                    Head[di] =!= DefiniteIntegrateZeroToInfinity,
                    di,
                    Integrate[
                      Collect[v Exp[-v^2/(4 t)](ff /. t->v),  Delta],
                      {v, 0, Infinity}
                    ] //. DeltaRules
                  ]
                ] / (2 Sqrt[Pi t^3])
              ] /;
              FreeQ[ff, InverseLaplace]
            ]
          ] /;
          (
            (SameQ @@ g) &&
            (
              MatchQ[g[[1]], _Times] ||
              MatchQ[g[[1]], _Plus]
            )
          ) || (And @@ (MatchQ[#, Power[s, Rational[_Integer, 2]]]& /@ g))
        ]
    },
    gg /;
    FreeQ[gg, InverseLaplace] && FreeQ[gg, Fail]
  ]


(************************************************************************)
End[]             (* end `Private` Context                              *)
(************************************************************************)
(************************************************************************)
EndPackage[]      (* end package Context                                *)
(************************************************************************)
(*:Limitations:

As in the case of integration, there are more special cases
which appear in the reference book, yet are not covered here.

One larger (and potentially needed) area which is not covered
is the treatment of functions which are defined differently
on different regions. Implementing those in MMa (e.g. using
If) is possible, but would require a wide support
(e.g. If[x<0, f1, 0] + If[x<0, 0, f2] -> If[x<0, f1, f2] etc.)
not available now.
Laplace can be evaluated numerically by performing the
integral in the definition using NIntegrate.
Although that is possible in principle for InverseLaplace,
the convergence of the typical integrals one encounters is
much poorer (goes like Sin[t]/t). If, at some point,
NIntegrate will be more successful in evaluating these,
a numeric rule for InverseLaplace could be added.

*)

(*:Tests:

Test[

Laplace[Sin[t], t, s],
(1 + s^2)^(-1)

]

Test[

InverseLaplace[%, s, t],

Sin[t]

]

NTest[

N[Laplace[1/(1+Exp[t]), t, 2.3]],

0.1735632309153487

]

*)
