(*  :Title:	Inverse Laplace Transform  *)

(*  :Authors:	Wally McClure, Brian Evans, James McClellan  *)

(*
    :Summary:	Multidimensional, bilateral inverse Laplace transform
		for signal processing expressions.
 *)

(*  :Context:	SignalProcessing`Analog`InvLaPlace`  *)

(*  :PackageVersion:  2.6	*)

(*
    :Copyright:	Copyright 1990-1991 by Brian L. Evans
		Georgia Tech Research Corporation

	Permission to use, copy, modify, and distribute this software
	and its documentation for any purpose and without fee is
	hereby granted, provided that the above copyright notice
	appear in all copies and that both that copyright notice and
	this permission notice appear in supporting documentation,
	and that the name of the Georgia Tech Research Corporation,
	Georgia Tech, or Georgia Institute of Technology not be used
	in advertising or publicity pertaining to distribution of the
	software without specific, written prior permission.  Georgia
	Tech makes no representations about the suitability of this
	software for any purpose.  It is provided "as is" without
	express or implied warranty.
 *)

(*
    :History:	date started -- April 25, 1990 (adapted from "InvZTransform.m")
		added Dialogue ability -- September, 1990
		handles real-valued polynomials -- March, 1991
 *)

(*  :Keywords:	Laplace transform, region of convergence	*)

(*  :Source:	Operational Mathematics by Churchill  *)

(*
    :Warning:	Works only on one dimensional and separable
		multidimensional transforms.
*)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*
    :Discussion:  1-D rules base has 70 rules:

		    I.   rational transforms       26 rules *
		    II.  non-rational transforms   28 rules
		    III. transform properties      10 rules *
		    IV.  transform SP structures    0 rules
		    V.   transform strategies       6 rules

		  * a partial fractions strategy rule exists in section

	Unlike some symbolic transform packages, these rules are maintained
	in  a  list  called  InvLaPlaceRules.   This  was necessary because
	Mathematica reordered  these  rules  if they were coded as a set of
	transformation functions.  Another benefit of the list form is that
	it allowed more control over how rules are applied.

	The  InvLaPlace  object (function) first  calls the interface rules
	(InvLaPlaceInterfaceRules), which  handle  default arguments, error
	messages, and multidimensional transforms.  These rules simply call
	MyInvLaPlace  once  per  dimension  of the transform.  This inverse
	transform is biased toward causal systems.   However, by specifying
	the proper  region of convergence (ROC),  anti-causal functions can
	be obtained (see InvLaPlace::usage).

	The  driver  for  the one-dimensional rule base is the six argument
	version  of  myinvlaplace.    myinvlaplace[X, s, t, rm, rp, st, op]
	either returns  the  inverse transform as a function of time (t) or
	an approximation to the actual inverse transform.  If the rule base
	can do neither,  which  should  not  happen, the interface function
	cleanup will report any errors.  Arguments of myinvlaplace:

	X  laplace-transform function		rm  Rminus of ROC
	s  laplace-transform variable		rp  Rplus of ROC
	t  "time"-domain variable		st  local state semaphores

	At  each  step  in the transform rule base,  the current expression
	has a local state  associated  with  it.   This state consists of a
	list of six boolean values.   Each boolean value is associated with
	a strategy.  If an element is True, then that strategy has not been
	tried yet; if False, then that strategy has already been tried, and
	it will not be tried again.   Thus,  local state is used to prevent
	infinite loops which would be caused by the  repetitive application
	of certain strategy rules.  See the section below called  S T A T E
	D E F I N I T I O N  and see section VI of the rules.
 *)

(*  :Functions:	 InvLaPlace  *)


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ];


(*  B E G I N     P A C K A G E  *)

BeginPackage [	"SignalProcessing`Analog`InvLaPlace`",
		"SignalProcessing`Analog`LSupport`",
		"SignalProcessing`Support`TransSupport`",
		"SignalProcessing`Support`SigProc`",
		"SignalProcessing`Support`ROC`",
		"SignalProcessing`Support`SupCode`",
		"SignalProcessing`Support`FilterSupport`" ]

(*  B E G I N     P A C K A G E  *)
		
(*  U S A G E     I N F O R M A T I O N  *)

InvLaPlace::usage =
	"InvLaPlace[f, s] and InvLaPlace[f, s, t] gives the multidimensional \
	bilateral inverse Laplace transform of f. \
	A region of convergence can be specified by using \
	InvLaPlace[{f, rm, rp}, s, t], where rm is R- and rp is R+ \
	in the region (strip) of convergence: R- < Re(s) < R+. \
	Note that InvLaPlaceTransform is an alias for InvLaPlace."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin[ "`Private`" ]


(*  U S E R     I N T E R F A C E  *)

(*  InvL operator  *)
Unprotect[InvL]
InvL/: TheFunction[ InvL[ s_, t_ ][ f_ ] ] := InvLaPlace[ f, s, t ]
Protect[InvL]

(*  InvLaPlace  *)
Options[InvLaPlace] :=
	{ Apart -> Rational, Definition -> False,
	  Dialogue -> True, Simplify -> True,
	  Terms -> 10, TransformLookup -> {} }

(*  bridge from outside world to interface to rule base  *)
InvLaPlace[ x_ ] :=
	invlaplacedriver[ Options[InvLaPlace], x ]
InvLaPlace[ x_, s_ ] :=
	invlaplacedriver[ Options[InvLaPlace], x, s ]
InvLaPlace[ x_, s_, t_ ] :=
	invlaplacedriver[ Options[InvLaPlace], x, s, t ]
InvLaPlace[ x_, s_, t_, op__ ] :=
	invlaplacedriver[ ToList[op] ~Join~ Options[InvLaPlace], x, s, t ]

(*  driver for interface to rule base  *)
invlaplacedriver[ op_, rest__ ] :=
	Block [	{},
		tlist = {};			(* global variable *)
		cleanup[ invlaplace[op, rest],
			 Replace[Dialogue, op],
			 TrueQ[Replace[Simplify, op]] ] ]

(*  interface to rule base   *)
invlaplace[ op_, args__ ] :=
	Block [ {},
		Replace[ linverse[op, args], InvLaPlaceInterfaceRules ] ]

InvLaPlaceInterfaceRules = {
	linverse[ op_, e_ ] :>
		invlaplace[ op, e, LVariables[e] ] /; LTransformQ[e],
	linverse[ op_, e_ ] :>
		Message[ InvLaPlace::novariables ],

	linverse[ op_, e_, s_ ] :>
		invlaplace[ op, e, s, DummyVariables[Length[s], Global`t] ] /;
		validvarQ[s],

	linverse[ op_, e_, s_, rest___ ] :>
		Message[ Transform::badvar, "s", InvLaPlace, s ] /;
		! validvarQ[s],
	linverse[ op_, e_, s_, t_, rest___ ] :>
		Message[ Transform::badvar, "time", InvLaPlace, t ] /;
		! validvarQ[t],

	linverse[ op_, e_List, s_, t_ ] :>
		invlaplace[ op, MakeLObject[e, s], s, t ],
	linverse[ op_, e_, s_Symbol, t_ ] :>
		MyInvLaPlace[ e, s, t, op ],
	linverse[ op_, e_, s_List, t_ ] :>
		multiDInvLaPlace[ e, s, t, op ]
}

cleanup[ trans_, dialogue_, simplify_ ] :=
	Block [	{adjtrans, lt, sdialogue},
		sdialogue = SameQ[dialogue, All];
		lt = If [ simplify,
			  SPSimplify[trans,
				     Dialogue -> sdialogue,
				     Variables -> tlist],
			  trans ];
		If [ InformUserQ[dialogue], Scan[explain, lt, Infinity] ];
		lt ]

explain[ invlaplace[ f_, s_, rest__ ] ]:=
	Message[ Transform::incomplete, "inverse Laplace transform", f, s ]

validvarQ[ x_Symbol ] := True
validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
validvarQ[ x_ ] := False


(*  M E S S A G E S  *)

InvLaPlace::badROC =
	"Improper region of convergence in ``."
InvLaPlace::badterms =
	"Option Terms is not an integer (Terms -> ``)."
InvLaPlace::novariables =
	"No Laplace variables specified in inverse Laplace transform."


(*  S T A T E     D E F I N I T I O N  *)

Lfactorfield = 1		(* factor denominator			  *)
Lexpandfield = 2		(* apply Expand[] to expression		  *)
Lexpandallfield = 3		(* apply ExpandAll[] to expression	  *)
Lthefunfield = 4		(* apply TheFunction to expression	  *)
Lapartfield = 5			(* apply Apart[] to expression		  *)
Lmyapartfield = 6		(* apply MyApart[] to rational poly.	  *)
Lnormalizefield = 7		(* normalize denominator		  *)
Lseriesfield = 8		(* approximate function with power series *)
Lsplitfield = 9			(* split rat. fun. into rat. poly & other *)
Ldefinitionfield = 10		(* apply defintion of inv. Laplace trans. *)

Lstatevariables = 10

initInvLstate[] := Table[True, {Lstatevariables}]
nullInvLstate[] := Table[False, {Lstatevariables}]


(*  S U P P O R T I N G     R O U T I N E S  *)

(*  Causal or anti-causal  *)
leftsided[ p_, t_ ] :=
	( p /. { CStep[ b_. t + c_. ] -> CStep[ -b t + c ], t -> - t } )

rightsided[ p_, t_ ] := p

leftOrRightSided[ a_, p_, t_, rp_ ] := rightsided[ p, t ] /; InfinityQ[rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ]  /; SameQ[Re[a], rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ]  /; SameQ[-Re[a], rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ]  /; N[Re[a] > rp]
leftOrRightSided[ a_, p_, t_, rp_ ] := rightsided[ p, t ]
 
(*  MyInvLaPlace  *)
MyInvLaPlace[ e_, s_, t_, op_ ] :=
	Block [ {rminus, rplus, valid},
		rminus = SPSimplify[ GetRMinus[e] ];
		rplus = SPSimplify[ GetRPlus[e] ];
		valid = Apply[And, Thread[ToList[rminus] < ToList[rplus]]];
		If [ TrueQ[! valid],
		     Message[ InvLaPlace::badROC, e ],
		     MyInvLaPlace[ TheFunction[e], s, t, rminus,
		     		   rplus, initInvLstate[], op ] ] ] /;
	LTransformQ[e]

MyInvLaPlace[ e_, s_, t_, op_ ] :=
	MyInvLaPlace[ e, s, t, -Infinity, Infinity, initInvLstate[], op ]

MyInvLaPlace[ e_, s_, t_, rm_, rp_, st_, op_ ] :=
	Block [	{laste, newe, newinvrules, norme, trace},

		(* add the time variable to the global list of time vars *)
		AppendTo[tlist, t];

		(* generate the rules *)
		newinvrules = TransformFixUp[ s, t, op, myinvlaplace, False,
					      InvLaPlace, Null, Null ];
		InvLaPlaceRules = Join[ newinvrules, OriginalInvLaPlaceRules ];

		(* determine the level of dialogue and pre-process f(t) *)
		trace = SameQ[ Replace[Dialogue, op], All ];
		norme = normalizeExpression[e, s];
		If [ trace && ! SameQ[norme, e],
		     Print[e]; Print["is first rewritten as"];
		     Print[norme]; Print["Then, its inverse transform is"] ];

		(* take the 1-D inverse Laplace transform *)
		newe = myinvlaplace[ norme, s, t, rm, rp, st, fixUp[op] ];
		While [ ! SameQ[laste, newe],
			If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
			laste = newe;
			newe = MapAll[transform, laste] ];

		(* fix up the result *)
		newe = postinvtransform[laste];
		While [ ! SameQ[laste, newe],
			If [ trace, Print[ newe]; Print[ "which becomes" ] ];
			laste = newe;
			newe = postinvtransform[laste] ];
		If [ trace, Print[ newe] ];

		(* Simplification rule based on t being real variable *)
		(* Mathematica 2.0 always assumes var to be complex   *)
		If [ TrueQ[ $VersionNumber >= 2.0 ],
		     newe = newe /.  ( ((t^l_.) x_)^k_ :> t^(l k) x^k /;
					( Positive[l] || Negative[l] ) &&
					( Positive[k] || Negative[k] ) ) ];

		newe /. ( f_ Delta[t + a_.] :> (f /. t -> -a) Delta[t + a] ) ]


(*  Kludge around the way Mathematica does Apart			*)
(*  If the user has selected Apart -> All, then use approximate roots	*)
(*  Otherwise, inverse transform it as a filter structure		*)
(*  Enable denominator normalization just in case Apart doesn't do it	*)

breakUpRealValuedPoly[x_, s_, t_, rm_, rp_, st_, op_] :=
	Block [	{dialogue, input, numer, state,},
		state = SetStateField[st, Lmyapartfield, False];
		state = SetStateField[state, Lnormalizefield, True];
		dialogue = InformUserQ[op];
		If [ dialogue,
		     Print[ "( since the Apart option is not All," ];
		     Print[ "  the inverse Laplace transform of" ];
		     Print[ "  ", x ];
		     Print[ "  will be an IIR filter whose input is" ];
		     Print[ "  the inverse transform of the numerator. )" ] ];

		numer = Numerator[x];
		input = If [ realValuedPolynomialQ[numer, s] &&
	  		       realValuedCoefficientsQ[numer, s],
			     CFIR[t, CoefficientList[numer, s]],
			     myinvlaplace[numer, s, t, rm, rp, state, op] ];
		CIIR[t, CoefficientList[Denominator[x], s]][input] ] /;
	! SameQ[ Replace[Apart, op], All ]

breakUpRealValuedPoly[x_, s_, t_, rm_, rp_, st_, op_] :=
	partialFractionsKludge[x, s, t, rm, rp, st, op, N, "approximate"]

(*  definitionDialogue  *)
definitionDialogue[ oldexpr_, trans_, operator_, options_ ] :=
	Block [	{},
		If [ InformUserQ[options],
		     Print[ "( after using ", operator,
			    "  to apply the definition to" ];
		     Print[ "  ", oldexpr ];
		     Print[ "  to find its transform"];
		     Print[ "  ", trans, " )" ] ];
		trans ]

(*  fixUp  *)
fixUp[ op_ ] := { Apart -> Replace[Apart, op],
		  Definition -> Replace[Definition, op],
		  Dialogue -> Replace[Dialogue, op],
		  Simplify -> Replace[Simplify, op],
		  Terms -> Replace[Terms, op] }


(*  multiDInvLaPlace  *)
multiDInvLaPlace[ e_, s_, t_, op_ ] :=
	MultiDInvTransform [ e, s, t, op, LTransformQ,
			     MyInvLaPlace, MakeLObject, Global`t ]

(*  normalizePoly and normalizeExpression  *)
truepoly[x_, s_] := (! FreeQ[x, s]) && PolynomialQ[x, s]

unnormalizedq[x_, s_] := ! SameQ[1, Coefficient[x, s, Exponent[x, s]]]

normalizePoly[poly_, s_, k_] :=
	Block [	{coeff, maxpower},
		maxpower = Exponent[poly, s];
		coeff = Coefficient[poly, s, maxpower];
		coeff^k ( Distribute[poly / coeff]^k ) ]

normalizeExpression[x_, s_] :=
	Block [ {},
		x /. { (poly_Plus)^k_. :> normalizePoly[Expand[poly], s, k] /;
			 truepoly[poly, s],
		       poly_ :> Expand[poly] /; truepoly[poly, s] } ]

(*  normalizeRatPoly  *)
normalizeRatPoly[x_, s_] := normalizeRatPoly[Numerator[x], Denominator[x], s]

normalizeRatPoly[numer_, denom_^k_., s_] :=
	Block [ {coeff, maxpower},
		maxpower = Exponent[denom, s];
		coeff = Coefficient[denom, s, maxpower];
		coeff^k numer / ( Distribute[denom / coeff]^k ) ]

(*  partialFractionsDialogue  *)
partialFractionsDialogue[p_, s_, options_] :=
	partialFractionsDialogue[ p, s, options,
				  normalizeExpression[Apart[p, s], s] ]

partialFractionsDialogue[p_, s_, options_, partfrac_, rootinfo_:"exact"] :=
	Block [ {dialogue},
		dialogue = InformUserQ[options];

		(*  Tell the user that we're applying partial fractions  *)
		If [ dialogue && ! SameQ[p, partfrac],
		     Print["( after performing partial fraction expansion on"];
		     Print["  ", p];
		     Print["  using its ", rootinfo, " roots:" ];
		     Print["  ", partfrac, " . )" ] ];

		partfrac ]

(*  partialFractionsKludge -- use MyApart because Apart has failed us  *)
partialFractionsKludge[ x_, s_, t_, rm_, rp_, st_, op_,
			filter_:Identity, roots_:"exact" ] :=
	Block [	{partfrac, result, state},
		state = SetStateField[st, Lapartfield, False];
		state = SetStateField[state, Lmyapartfield, False];
		state = SetStateField[state, Lnormalizefield, False];
		partfrac = normalizeExpression[MyApart[x, s, filter], s];
		result = partialFractionsDialogue[x, s, op, partfrac, roots];
		myinvlaplace[ result, s, t, rm, rp, state, op ] ]

(*  postinvtransform  *)
postinvtransform[f_] := ( f /. postinvrules )

postinvrules = {
	derivative[f_, t_Symbol, k_, False] :>
		D[ postinvtransform[f], {t, k} ],
	derivative[f_, t_Symbol, k_, v_] :>
		Block [	{curderivative, ffun, i, limit, result},
			ffun = postinvtransform[f];
			result = Unit[k-1][t] Limit[ffun, t -> v];
			curderivative = ffun;
			For [ i = 2, i <= k, i++,
			      curderivative = D[curderivative, t];
			      limit = Limit[curderivative, t -> v];
			      result += Unit[k - i][t] limit ];
			result + D[curderivative, t] ],
	integrate[a_ + b_, t_, dummy_Symbol] :>
		integrate[postinvtransform[a], t, dummy] +
		integrate[postinvtransform[b], t, dummy],
	integrate[f_. CStep[dummy_ + b_.], t_, dummy_Symbol ] :>
		Integrate[f, {dummy, -b, t}] CStep[t + b],
	integrate[f_. CStep[a_. dummy_ + b_.], t_, dummy_Symbol ] :>
		( CStep[a t + b] / a ) *
		Integrate[f /. t -> t / a,
			  {dummy, Max[-b, -Infinity a], a t}],
	substitutefort[f_, t_, newt_] :>
		( postinvtransform[f] /. t -> newt ),

	(*  These two rules are needed to support a series 	*)
	(*  expansion which returns a constant term plus	*)
	(*  terms of the form  approx * s^r   .			*)
	myinvlaplace[ a_, s_, t_, rm_, rp_, st_, op_] :> a Delta[t] /;
		FreeQ[a, s],		   (* Converges all ROC *)

	myinvlaplace[ a_. s_^r_., s_, t_, rm_, rp_, st_, op_] :>
		a Unit[r][t] /;
		FreeQ[{a,r}, s],	   (* Converges all ROC *)

	(*  Attempt a series expansion about s = 0		*)
	(*  This is the strategy when all else has failed.	*)
	(*  Please keep this as the last rule.			*)
	myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
		SeriesStrategy[ x, s, t, rm, rp, st, op ] /;
		GetStateField[st, Lseriesfield],

	x_ :> x
}

(*  realAndPosPart  *)
realposcheckq[x_] := Implies[ NumberQ[x], RealQ[x] && x > 0 ]
realnegcheckq[x_] := Implies[ NumberQ[x], RealQ[x] && x < 0 ]

realAndPosPart[0] := 0
realAndPosPart[x_Times] := Map[realAndPosPart, x]
realAndPosPart[Complex[a_Integer, b_Integer]] := realAndPosPart[GCD[a,b]]
realAndPosPart[Complex[b_, b_]] := realAndPosPart[b]
realAndPosPart[Complex[0, b_]] := realAndPosPart[b]
realAndPosPart[x_Symbol] := x
realAndPosPart[x_] := x /; realposcheckq[ N[x] ]
realAndPosPart[x_] := -x /; realnegcheckq[ N[x] ]
realAndPosPart[x_] := 1

(*  realValuedCoefficientsQ  *)
realValuedCoefficientsQ[z_, t_] :=
	Apply[Or, Map[RealQ, CoefficientList[z, t]]]

(*  realValuedPolynomialQ  *)
realValuedPolynomialQ[z_] := realValuedPolynomialQ[z, Global`x]

realValuedPolynomialQ[c_. t_^r_., t_] := IntegerQ[r] && r > 0 && RealValuedQ[c]
realValuedPolynomialQ[c_, t_] := RealValuedQ[c]
realValuedPolynomialQ[x_ + y_, t_] :=
	realValuedPolynomialQ[x, t] && realValuedPolynomialQ[y, t]

(*  seriesDialogue  *)
seriesDialogue[x_, expansion_, terms_, options_] :=
	Block [	{dialogue},
		dialogue = InformUserQ[options];

		If [ dialogue,
		     Print["( after breaking up the expression"];
		     Print["  ", x];
		     Print["  into its series representation"];
		     Print["  ", expansion];
		     Print["  the inverse Laplace transform will now be",
			   "  applied to the first ", terms, " terms. )"] ] ]

(*  SeriesExpansion  *)
SeriesExpansion[ True, x_, s_, power_ ] := Series[x, {s, 0, power}]

SeriesExpansion[ False, x_, s_, power_ ] :=
	Block [	{negx, expansion},
		negx = x /. s -> s^-1;
		expansion = Series[negx, {s, 0, power}];
		expansion /. s -> s^-1 ]

(*  SeriesStrategy  *)
SeriesStrategy[ x_, s_, t_, rm_, rp_, st_, op_ ] :=
	Block [ {expansion, exponents, posexpandflag, state, terms},
		terms = Replace[Terms, op];
		state = SetStateField[st, Lseriesfield, False];

		(*  If the Terms option is disabled, then take no action.  *)
		If [ ! terms,
		     Return[ myinvlaplace[x, s, t, rm, rp, state, op] ] ];

		If [ ! IntegerQ[terms],
		     Message[ InvLaPlace::badterms, terms ];
		     Return[ myinvlaplace[x, s, t, rm, rp, state, op] ] ];

		(*  See if we should expand in pos. or neg. powers	*)
		exponents = GetAllExponents[x, s];
		posexpandflag = If [ Apply[And, Thread[exponents >= 0]],
				     True, False, False ];

		(*  First expansion  *)
		expansion = SeriesExpansion[posexpandflag, x, s, terms - 1];

		(*  Second expansion if first failed  *)
		If [ ! SameQ[Head[expansion], SeriesData],
		     expansion = SeriesExpansion[Not[posexpandflag], x, s, terms - 1] ];

		(*  If this fails, call MyInvLaPlace again		*)
		(*  but with series expansion disabled			*)
		If [ ! SameQ[Head[expansion], SeriesData],
		     myinvlaplace[ x, s, t, rm, rp, state, op ],
		     seriesDialogue[x, expansion, terms, op];
		     state = nullInvLstate[];
		     Map [ myinvlaplace[#1, s, t, rm, rp, state, op]&, 
			   Normal[expansion] ] ] ]

(*  similarityQ  *)
similarityScale[f_, s_] := realAndPosPart[ ScalingFactor[f, s] ]
similarityQ[f_, s_] :=
	Block [	{scale},
		scale = similarityScale[f, s];
		(! SameQ[scale, 1]) && (! SameQ[scale, 0]) ]

(*  splitratfun *)
keepterm[ x_ + y_, s_ ] := keepterm[x, s] + keepterm[y, s]
keepterm[ a_. s_^n_., s_ ] := a s^n /; FreeQ[{a, n}, s]
keepterm[ a_, s_ ] := a /; FreeQ[a, s]
keepterm[ a_, s_ ] := 0

splitratfun[x_, s_] := splitratfun[ Numerator[x], Denominator[x], s ]

splitratfun[ numer_, denom_, s_ ] :=
	Block [ {expnumer, polyexpr, restexpr},
		expnumer = Expand[numer];
		polyexpr = keepterm[expnumer, s];
		restexpr = expnumer - polyexpr;
		{ polyexpr / denom, restexpr / denom } ]

(*  transform  *)
transform[ expr_ ] :=
	If [ SameQ[Head[expr], myinvlaplace],
	     Replace[expr, InvLaPlaceRules],
	     expr ]


Format[ myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] ] :=
	SequenceForm[ ColumnForm[{"L" Superscript[-1],
				  "  " ~StringJoin~ ToString[s]}],
		      { x } ]

Format[ derivative[x_, t_Symbol, k_, flag_] ] := 
	SequenceForm[ ColumnForm[{"D" Superscript[k],
				  "  " ~StringJoin~ ToString[t]}],
		      { x } ]

Format[	substitutefort[f_, t_, newt_] ] :=
	SequenceForm[ {f}, Subscript[t -> newt] ]

Format[ integrate[f_. CStep[a_. dummy_ + b_.], t_, dummy_Symbol ] ] :=
	SequenceForm[ "Integrate[", f, ", ", {dummy, 0, t}, "] ",
		      CStep[a t + b] ]

Format[ integrate[f_, t_, dummy_Symbol ] ] :=
	SequenceForm[ "Integrate[", f, ", ", {dummy, 0, t}, "]" ]


(*  B E G I N     R U L E S  *)


InvLaPlaceRules = {}


OriginalInvLaPlaceRules = { 


  (* I.   Rational Transform Pairs				*)

  (*	  A.  constant functions				*)
  myinvlaplace[ a_, s_, t_, rm_, rp_, st_, op_] :> a Delta[t] /;
	FreeQ[a, s],					(* Converges all ROC *)


  (*	  B.  constant times exponential			*)
  myinvlaplace[ a_. Exp[b_. s_], s_, t_, rm_, rp_, st_, op_ ] :>
	a Delta[t + b] /;
	FreeQ[{a, b}, s],
 

  (*	  C.  inverse of sinusoidal forms			*)
  myinvlaplace[ c_. ( w_. Cos[w_] + s_ Sin[w_] ) / (s_^2 + w_^2), s_, t_,
		rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ 0, c Sin[w + w t] CStep[t], t, rp ] /;
	FreeQ[{c, w}, s],

  myinvlaplace[ c_. ( s_ Cos[b_] - w_. Sin[b_] ) / ( s_^2 + w_^2 ), s_, t_,
		rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ 0, c Cos[b + w t] CStep[t], t, rp ] /;
	FreeQ[{w,b,c}, s],


  (*	  D.  sections that are functions of single power of s	*)
  (*	      1.  single pole sections				*)
  myinvlaplace[ b_. (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ a,
			  b t^-(n+1) Exp[- a t] CStep[t] / Gamma[-n],
			  t, rp ] /;
	FreeQ[{a,b,n}, s] && Negative[n],
 
  myinvlaplace[ b_. / (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ a,
			  b t^(n-1) Exp[- a t] CStep[t] / Gamma[n],
			  t, rp ] /;
	FreeQ[{a,b,n}, s],

  (*	      2.  single zero sections				*)
  myinvlaplace[ b_. (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
	b Exp[- a t] Unit[n][t] /;
	FreeQ[{a,b,n}, s],



  (*	  D.  second-order sections				*)
  (*          only encode the cosine/sine pairs because Mma	*)
  (*          will automatically convert Sin[I x] to I Sinh[x]	*)
  (*	      and Cos[I x] to Cosh[x].				*)
  myinvlaplace[ b_. / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[Sqrt[a], b Sin[Sqrt[a] t] CStep[t] / Sqrt[a], t, rp] /;
	FreeQ[{a,b}, s],

  myinvlaplace[ b_. s_ / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ Sqrt[a], b Cos[Sqrt[a] t] CStep[t], t, rp ] /; 
	FreeQ[{a,b}, s],

  myinvlaplace[ (s_ + a_)/(s_^2 + c_. s_ + d_), s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{b, invtrans},
		b = Sqrt[d - a^2];			(* d = a^2 + b^2 *)
		invtrans = Exp[-a t] Cos[b t] CStep[t];
		leftOrRightSided[-a, invtrans, t, rp] ] /;
	FreeQ[{a,c,d}, s] && SameQ[c, 2 a],

  myinvlaplace[ k_. / ( s_^2 + c_. s_ + d_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{a},
		a = c/2;		(* complete the square *)
		Exp[-a t] myinvlaplace[ k / ( s^2 + d - a^2 ), s, t,
					rm, rp, st, op ] ] /;
	FreeQ[{c,d,k}, s],


  (*	  E.  third order sections				*)
  myinvlaplace[ b_. /( s_^3 + a_), s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ Max[-Re[(-a)^(1/3)/2], -Re[a^(1/3)]],
			  b / (3 a^(2/3)) (Exp[- a^(1/3) t] - Exp[a^(1/3) t/2]
				( Cos[ a^(1/3) t Sqrt[3] / 2 ] - 
				  Sqrt[3] Sin[ a^(1/3) t Sqrt[3] / 2 ] ) )
			  CStep[t],
			  t,
			  rp ] /;
	FreeQ[{a, b}, s],				(* [Churchill] *)


  (*	  F.  fourth order sections				*)
  myinvlaplace[ b_. / ( s_^2 + a_ )^2, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{k},
		k = Sqrt[a] t;
		leftOrRightSided[ Sqrt[-a],
			   	  b (Sin[k] - k Cos[k]) CStep[t] / (2 a^(3/2)),
				  t, rp ] ] /;
	FreeQ[{a, b}, s],

  myinvlaplace[ b_. s_ / ( s_^2 + a_ )^2, s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ Sqrt[-a],
			  b t Sin[Sqrt[a] t] CStep[t] / (2 Sqrt[a]),
			  t, rp ] /;
	FreeQ[{a, b}, s],

  myinvlaplace[ b_. / ( s_^4 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{k},
		k = (a/4)^(1/4);
		leftOrRightSided[ Re[ a^(1/4) / Sqrt[2] ],
				  b (Sin[k t] Cosh[k t] - Cos[k t] Sinh[k t])
				    CStep[t] / ( 4 k^3 ),
				  t,
				  rp ] ] /;
	FreeQ[{a, b}, s],

  myinvlaplace[ b_. s_ / ( s_^4 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{k},
		k = (a/4)^(1/4);
		leftOrRightSided[ Re[ a^(1/4) / Sqrt[2] ], 
				  b Sin[k t] Sinh[k t] CStep[t] / ( 2 k^2 ),
				  t,
				  rp ] ] /;
	 FreeQ[{a, b}, s],


  (*	  G.  higher order sections				*)
  myinvlaplace[ b_. s_^n_., s_, t_, rm_, rp_, st_, op_ ] :>
  	leftOrRightSided[ 0, b t^-(n + 1) / Gamma[-n] CStep[t], t, rp ] /; 
	Negative[n] && FreeQ[b, s],

  myinvlaplace[ b_. s_^n_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{index = Unique["index"]}, 
		leftOrRightSided[ 0,
			   b CStep[t] 2^(-n-1/2) t^(-n-1) /
			   ( Product[index, {index, 1, -2 n - 1}] Sqrt[Pi] ),
			   t, rp ] ] /; 
	IntegerQ[ - n - 1/2 ] && FreeQ[b, s],

  myinvlaplace[ b_. ( s_^2 + a_ )^k_., s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ a,
			  b Sqrt[Pi] t^(-k - 1/2) / ((2 Sqrt[a])^(-k - 1/2))
			    BesselJ[-k - 1/2, Sqrt[a] t ] CStep[t] / Gamma[-k],
			  t, rp ] /;
	FreeQ[{a, b}, s] && Negative[k],

  myinvlaplace[ b_. / ( s_^2 + a_ )^k_., s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ a,
			  b Sqrt[Pi] t^(k -1/2) / ((2 Sqrt[a])^(k - 1/2))
			    BesselJ[k - 1/2, Sqrt[a] t ] CStep[t] / Gamma[k],
			  t, rp ] /;
	FreeQ[{a,b,k}, s],


  (*	  H.  partial fractions of rational polynomials		*)

  (*          1. split expression into rational polys and rest  *)
  (*             this must occur before partial fractions	*)
  (*		 because Apart ONLY handles ratio of polys	*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{splitlist, state},
		splitlist = splitratfun[x, s];
		state = SetStateField[st, Lsplitfield, False];
		myinvlaplace[ First[splitlist], s, t, rm, rp, state, op] +
		  myinvlaplace[ Second[splitlist], s, t, rm, rp, state, op] ] /;
	GetStateField[st, Lsplitfield] && GetStateField[st, Lapartfield] &&
	  truepoly[Denominator[x], s] && (! PolynomialQ[Numerator[x], s]),

  (*          2. Reduce all rational polynomials into		*)
  (*	         a sum of single pole, second order sections,	*)
  (*	         and sections like 1 / ( s^4 + a^4 )		*)
  (*	         Enable denominator normalization just in case	*)
  (*	         partial fractions does not do it automatically	*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{state},
		state = SetStateField[st, Lapartfield, False];
		state = SetStateField[state, Lnormalizefield, True];
		myinvlaplace[ partialFractionsDialogue[x, s, op], s, t,
			      rm, rp, state, op ] ] /;
	GetStateField[st, Lapartfield] && RationalPolynomialQ[x, s],
 

  (*      I.  partial fractions that Apart could not resolve	*)
  (*	      perform MyApart to redo partial fractions for	*)
  (*	      denominators with real-valued coefficients	*)
  (*	      that are not handled by Apart. 		*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
	breakUpRealValuedPoly[x, s, t, rm, rp, st, op] /;
	GetStateField[st, Lmyapartfield] &&
	  realValuedPolynomialQ[Denominator[x], s] &&
	  realValuedCoefficientsQ[Denominator[x], s],




  (* II.  Non-rational Transform Pairs				*)

  (*	  A.  square root forms					*)
  myinvlaplace[ 1 / Sqrt[s_ + b_.], s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ 0, Exp[- b t] CStep[t] / (Sqrt[t] Sqrt[Pi]),
			  t, Infinity ] /;
	FreeQ[b, s],

  myinvlaplace[ s_ / ( s_ + a_ )^(3/2), s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ -Re[a], Exp[- a t] (1 - 2 a t) CStep[t] / Sqrt[Pi],
			  t, rp ] /;
	FreeQ[a, s],

  myinvlaplace[ 1 / ((s_ + a_.)^(3/2) Sqrt[s_ + b_.]), s_, t_,
		rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ a,
		t Exp[- (a + b) t / 2] ( BesselI[0, (a - b) t / 2] -
					 BesselI[1, (a - b) t / 2] ) CStep[t],
		t, rp ] /;
	FreeQ[{a,b}, s] && ! SameQ[a, b],

  myinvlaplace[ 1 / ( Sqrt[s_] + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	(1 / (Sqrt[Pi] Sqrt[t]) - a Exp[a^2 t](1 - Erf[a Sqrt[t]])) CStep[t] /; 
	FreeQ[a, s],

  (*	      general form of #38 and #39 [Churchill, 326],	*)
  (*	      Sqrt[s] / ( s + a ), augmented by s -> s + b	*)
  myinvlaplace[ Sqrt[s_ + b_.] / ( s_ + c_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  	leftOrRightSided[c,
			 ( 1 / (Sqrt[Pi] Sqrt[t]) - Sqrt[b - c] Exp[(b - c) t]
			   Erf[ Sqrt[(b - c) t] ] ) CStep[t],
			 t, rp ] /;
	FreeQ[{b,c}, s],

  (*	      general form of #40 and #41 [Churchill, 326],	*)
  (*	      1 / ( Sqrt[s] (s + a) ), augmented by s -> s + b	*)
  myinvlaplace[ 1 / (Sqrt[s_ + b_.] (s_ + c_)), s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[c,
			 Exp[- c t] Erf[ Sqrt[(b - c) t] ] CStep[t] /
			   Sqrt[b - c],
			 t, rp ] /;
	FreeQ[{b,c}, s],

  (*	      forms that inverse transform to Bessel functions	*)
  (*	      1/Sqrt[s^2 + a] covered as a higher order section	*)
  myinvlaplace[ ( Sqrt[s_^2 + a_] - s_ )^k_, s_, t_, rm_, rp_, st_, op_ ] :>
	leftOrRightSided[ a, k Sqrt[a]^k BesselJ[k, Sqrt[a] t] CStep[t] / t,
			  t, rp ] /;
	FreeQ[{a,k}, s] && Positive[k],

  myinvlaplace[ ( Sqrt[s_^2 + a_] - s_ )^v_. / Sqrt[s_^2 + a_], s_, t_,
		rm_, rp_, st_, op_ ] :>
	Block [	{timefun},
		Assuming[ Positive[v + 1], op ];
		timefun = Sqrt[a]^v BesselJ[v, Sqrt[a] t] CStep[t];
		leftOrRightSided[a, timefun, t, rp] ] /;
	FreeQ[{a,v}, s] && Implies[NumberQ[N[v]], N[v > -1]],

  myinvlaplace[ 1 / ( Sqrt[s_ + a_.] Sqrt[s_ + b_.] ), s_, t_,
		rm_, rp_, st_, op_ ] :>
	Exp[- (a + b) t / 2] BesselI[0, (a - b) t / 2] CStep[t] /;
	FreeQ[{a, b}, s],

  myinvlaplace[ ( s_ - Sqrt[s_^2 + a_] )^v_. / Sqrt[s_^2 + a_], s_, t_,
		rm_, rp_, st_, op_ ] :>
	Block [	{timefun},
		Assuming[ Positive[v + 1], op];
		timefun = Sqrt[-a]^v BesselI[v, Sqrt[-a] t] CStep[t];
		leftOrRightSided[ a, timefun, t, rp ] ] /;
	FreeQ[{a,v}, s] && Implies[NumberQ[N[v]], N[v > -1]],


  (*	  B.  hyperbolic forms 					*)
  myinvlaplace[ Tanh[k_. s_] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{n, context}, 
		context = $Context;
		$Context = "Global`";
		n = Unique["n"];			(* n is an integer *)
		$Context = context;
		(-1)^(n-1) ( CStep[t - 2 k (n - 1)] - CStep[t - 2 k n] ) ] /;
	FreeQ[k, s],

  myinvlaplace[ Coth[b_. s_] / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	Abs[Sin[Sqrt[a] t] ] / Sqrt[a] CStep[t] /;
	b == ( Pi / ( 2 Sqrt[a] ) ),


  (*	  C.  exponential forms					*)
  myinvlaplace[ Exp[a_. Sqrt[s_]] / Sqrt[s_], s_, t_,
		rm_, rp_, st_, op_ ] :>
	Exp[-a^2 / (4 t)] CStep[t] / (Sqrt[Pi] Sqrt[t]) /;
	FreeQ[a, s],

  myinvlaplace[ Exp[k_. / s_] s_^u_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [ {arg, defarg, exp, order},
		defarg = 2 Sqrt[k t];
		arg = 2 Sqrt[k] Sqrt[t];
		order = -u - 1;
		exp = order/2;
		Which [ Negative[k],
			  t^exp BesselJ[order, arg] CStep[t] / k^exp,
			Positive[k],
			  t^exp BesselI[order, arg] CStep[t] / k^exp,
			True,
			  t^exp BesselI[order, defarg] CStep[t] / k^exp ] ] /;
	Negative[u],

  myinvlaplace[ Exp[k_. / s_], s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{arg, scale},
		arg = 2 Sqrt[-k] Sqrt[t];
		scale = Sqrt[-k] / ( 2 Sqrt[t] );
		Delta[t] + scale ( BesselJ[-1, arg] - BesselJ[1, arg] ) ] /;
	FreeQ[k, s],

  myinvlaplace[ Exp[ k_ Sqrt[s_] ], s_, t_, rm_, rp_, st_, op_ ] :>
	- k Exp[- k^2 / (4 t)] CStep[t] / (2 Sqrt[Pi] Sqrt[t^3]) /;
	Negative[k],


  (*	  D.  logarithmic forms					*)
  myinvlaplace[ Log[c_. s_ + b_.] (c_. s_ + b_.)^k_., s_, t_,
		rm_, rp_, st_, op_ ] :>
   	Exp[- b t / Abs[c]] (t / Abs[c])^(-k - 1)
	  ( Gamma[-k] PolyGamma[-k] / Abs[Gamma[-k]]^2 -
	    Log[t / Abs[c]] / Gamma[-k] ) CStep[t] / Abs[c] /;
	FreeQ[{b,c}, s] && Negative[k],

  myinvlaplace[ Log[c_. s_ + b_.] / (c_. s_ + b_.)^k_., s_, t_,
		rm_, rp_, st_, op_ ] :>
   	Exp[- b t / Abs[c]] (t / Abs[c])^(k - 1)
	  ( Gamma[k] PolyGamma[k] / Abs[Gamma[k]]^2 -
	    Log[t / Abs[c]] / Gamma[k] ) CStep[t] / Abs[c] /;
	FreeQ[{b,c,k}, s],

  myinvlaplace[ Log[b_. s_ + a_] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
	( Log[a] - ExpIntegralEi[- a t / Abs[b]] ) CStep[t] /;
	FreeQ[{a,b}, s],

  myinvlaplace[ Log[b_. s_ + a_.], s_, t_, rm_, rp_, st_, op_ ] :>
  	- Abs[b] Exp[- a t / Abs[b]] CStep[t] / t /;
	FreeQ[{a,b}, s],

  myinvlaplace[ Log[a_. + b_. s_^2 ] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
	2 ( Log[Sqrt[a]] - CosIntegral[Sqrt[a] t / Sqrt[b]] ) CStep[t] /;
	 FreeQ[{a,b}, s],

  myinvlaplace[ Log[b_. s_] / (c_. s_^2 + 1), s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{time = t / Abs[b]},
		( Cos[time] SinIntegral[time] -
		  Sin[time] CosIntegral[time] ) CStep[t] / Abs[b] ] /;
	FreeQ[{b,c}, s] && c == b^2,

  myinvlaplace[ s_ Log[b_. s_] / (c_. s_^2 + 1), s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{time = t / Abs[b]},
		( - Cos[time] CosIntegral[time] -
		    Sin[time] SinIntegral[time] ) CStep[t] / (b Abs[b]) ] /;
	FreeQ[{b,c}, s] && c == b^2,

  myinvlaplace[ Log[ a_ / b_ ], s_, t_, rm_, rp_, st_, op_ ] :>
  	myinvlaplace[ Log[a] - Log[b], s, t, rm, rp, st, op ],

  myinvlaplace[ Log[ a_ b_ ], s_, t_, rm_, rp_, st_, op_ ] :>
  	myinvlaplace[ Log[a] + Log[b], s, t, rm, rp, st, op ],


  (*	  E.  arc tangent					*)
  myinvlaplace[ ArcTan[ k_ / s_ ], s_, t_, rm_, rp_, st_, op_ ] :>
  	Sin[ k t ] CStep[t] / t /; FreeQ[k,s],


  (*	  F.  BesselK						*)
  myinvlaplace[ BesselK[ 0, k_. s_ ], s_, t_, rm_, rp_, st_, op_ ] :>
   	1 / Sqrt[t^2 - k^2] CStep[t - k] /;
	FreeQ[k, s],

  myinvlaplace[ BesselK[-1, v_. Sqrt[s_]] / Sqrt[s_], s_, t_,
		rm_, rp_, st_, op_ ] :>
	Exp[-v^2/(4 t)] CStep[t] / v /;
	FreeQ[v, s],


  (* III. Transform Properties					*)


  (*	  A.  homogeneity -- pick off constants			*)
  myinvlaplace[ a_ f_, s_, t_, rm_, rp_, st_, op_ ] :> 
  	a myinvlaplace[ f, s, t, rm, rp, st, op ] /;
	FreeQ[a, s],

  (*	  -.  perform MyApart to redo partial fractions for	*)
  (*	      denominators with repeated roots that were not	*)
  (*	      properly handled by Apart.  kludge 		*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
	partialFractionsKludge[x, s, t, rm, rp, st, op] /;
	GetStateField[st, Lmyapartfield] && RationalPolynomialQ[x, s],


  (*	  D.  additivity					*)
  myinvlaplace[ x_ + y_, s_, t_, rm_, rp_, st_, op_ ] :>
  	myinvlaplace[ x, s, t, rm, rp, st, op ] + 
	myinvlaplace[ y, s, t, rm, rp, st, op ], 

  myinvlaplace[ (x_ + y_) / c_, s_, t_, rm_, rp_, st_, op_ ] :>
  	myinvlaplace[ x / c, s, t, rm, rp, st, op ] + 
	myinvlaplace[ y / c, s, t, rm, rp, st, op ], 

  myinvlaplace[ a_. Summation[i_, ib_, ie_, inc_][x_], s_, t_,
		rm_, rp_, st_, op_ ] :>
	Summation[i, ib, ie, inc][ myinvlaplace[ a x, s, t, rm, rp, st, op ] ],


  (*	  E.  multiplication by an exponential			*)
  myinvlaplace[ a_. Exp[b_. s_ + c_.], s_, t_, rm_, rp_, st_, op_ ] :>
	substitutefort[ myinvlaplace[ a Exp[c], s, t, rm, rp, st, op ],
			t, t + b ] /;
	FreeQ[b, s],


  (*	  F.  similarity					*)
  (*          only use real, positive scaling factors		*)
  myinvlaplace[ f_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [ {newf, newrm, newrp, scale},
		scale = similarityScale[f, s];
		Assuming[ Positive[scale], op];
		newf = f /. s -> s / scale;
		newrm = If [ SameQ[rm, -Infinity], rm, scale rm ];
		newrp = If [ SameQ[rp, Infinity], rp, scale rp ];
		substitutefort[ myinvlaplace[newf, s, t, newrm, newrp, st, op],
				t, t / scale ] / scale ] /;
	similarityQ[f, s],


  (*      G.  shift						*)
  myinvlaplace[ f_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{newrm, newrp, normexpr, shift},
		{shift, normexpr} = GetShiftFactor[f, s];
		normexpr = normexpr /. s -> (s - shift);
		newrm = rm + Re[shift];
		newrp = rp + Re[shift];
		Exp[- shift t] *
		  myinvlaplace[normexpr, s, t, newrm, newrp, st, op] ] /;
	! SameQ[ First[GetShiftFactor[f, s]], 0 ],


  (*	  H.  pole in denominator				*)
  myinvlaplace[ f_ / ( s_ + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	myinvlaplace[ partialFractionsDialogue[f / (s + a), s, op], s, t,
		      rm, rp, SetStateField[st, Lapartfield, False], op ] /;
	GetStateField[st, Lapartfield] && FreeQ[a, s],

  myinvlaplace[ f_ / ( s_ + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
	Exp[- a t] myinvlaplace[(f /. s -> s - a) / s, s, t, rm, rp, st, op] /;
	FreeQ[a, s],


  (*	  I.  multiplication by s is progressive		*)
  (*	      differentiation in the time domain;		*)
  (*	      progressiveness handled by postinvrules		*)
  myinvlaplace[ s_^k_. a_., s_, t_, -Infinity, Infinity, st_, op_ ] :>
	Block [ {dialogue},
                dialogue = SameQ[ Replace[Dialogue, op], All ];
		Assuming[ k > 0, dialogue ];
		If [ dialogue && ! IntegerQ[k],
		     Print[ "where ", k, " is an integer" ] ];
  		derivative[ myinvlaplace[a, s, t, rm, rp, st, op],
			    t, k, False ] ] /;
	FreeQ[k, s] && Implies[ NumberQ[k], IntegerQ[k] && ( k > 0 )],

  myinvlaplace[ s_^k_. a_., s_, t_, 0, Infinity, st_, op_ ] :>
	Block [ {dialogue},
                dialogue = SameQ[ Replace[Dialogue, op], All ];
		Assuming[ k > 0, dialogue ];
		If [ dialogue && ! IntegerQ[k],
		     Print[ "where ", k, " is an integer" ] ];
  		derivative[myinvlaplace[a, s, t, rm, rp, st, op], t, k, 0] ] /;
	FreeQ[k, s] && Implies[ NumberQ[k], IntegerQ[k] && ( k > 0 )],


  (*	  J.  division by s is integration in the time domain	*)
  myinvlaplace[ f_ s_^n_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{tau = Unique["tau"]},
		integrate[ myinvlaplace[f s^(n + 1), s, tau, rm, rp, st, op],
			   t, tau ] ] /;
	IntegerQ[n] && (n < 0),




  (* V.   Strategies						*)

  (*	  A.  partial fractions					*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
	myinvlaplace[ partialFractionsDialogue[x, s, op], s, t, rm, rp,
		      SetStateField[st, Lapartfield, False], op ] /;
	GetStateField[st, Lapartfield],


  (*      B. Normalize denominator polynomials			*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
	myinvlaplace[ normalizeRatPoly[x, s], s, t, rm, rp,
		      SetStateField[st, Lnormalizefield, False], op ] /;
	GetStateField[st, Lnormalizefield] && unnormalizedq[Denominator[x], s],


  (*	  C.  factor denominator				*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
   	Block [	{denom, numer},
		denom = Denominator[x];
		numer = Numerator[x];
		denom = Factor[denom];
		myinvlaplace[ numer / denom, s, t, rm, rp,
			      SetStateField[st, Lfactorfield, False], op ] ] /;
	GetStateField[st, Lfactorfield],


  (*	  D.  distribute expression				*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  	myinvlaplace[ Distribute[x], s, t, rm, rp, st, op ] /;
	SameQ[Head[x], Times] && ! SameQ[x, Distribute[x]],


  (*	  E.  expand numerators of expression			*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  	myinvlaplace[ Expand[x], s, t, rm, rp,
		      SetStateField[st, Lexpandfield, False], op ] /;
	GetStateField[st, Lexpandfield],


  (*	  F.  expand all terms in expression			*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  	myinvlaplace[ ExpandAll[x], s, t, rm, rp,
		      SetStateField[st, Lexpandallfield, False], op ] /;
	GetStateField[st, Lexpandallfield],


  (*	  G.  replace new operators and functions		*)
  (*	      with their mathematical definition		*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
       	myinvlaplace[ TheFunction[x], s, t, rm, rp,
		      SetStateField[st, Lthefunfield, False], op ] /;
	GetStateField[st, Lthefunfield],

 
  (*      H.  apply definition if the user has enabled that	*)
  (*	      the Definition option.				*)
  myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
	Block [	{state, trans},
		state = SetStateField[st, Ldefinitionfield, False];
		trans = Integrate[x Exp[s t], {s, rm, rp} ];
		If [ SameQ[Head[trans], Integrate],
		     myinvlaplace[x, s, t, rm, rp, state, op],
		     definitionDialogue[x, trans/(2 Pi), Integrate, op] ] ] /;
	GetStateField[st, Ldefinitionfield] && Replace[Definition, op]

}

(*  E N D     R U L E S  *)


(*  E N D     P A C K A G E  *)

End[]
EndPackage[]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ];


(*  A L I A S E S  *)

Unprotect[ { InverseLaPlace, InverseLaPlaceTransform } ]
InverseLaPlace = InvLaPlace
InverseLaPlaceTransform = InvLaPlace
Protect[ { InverseLaPlace, InverseLaPlaceTransform } ]


(*  H E L P     I N F O R M A T I O N  *)

Combine[ SPfunctions, { InvLaPlace } ]
Protect[ InvLaPlace ]


(*  E N D I N G     M E S S A G E  *)

Print[ "The inverse Laplace transform rule base InvLaPlace has been loaded." ]
Null
