(*^

::[paletteColors = 128; 
	fontset = title, "New York", 24, L3, center, bold, nohscroll;
	fontset = subtitle, "New York", 18, L2, center, bold, nohscroll;
	fontset = subsubtitle, "New York", 14, L2, center, bold, nohscroll;
	fontset = section, "New York", 14, L2, bold, nohscroll, grayBox;
	fontset = subsection, "New York", 12, L2, bold, nohscroll, blackBox;
	fontset = subsubsection, "New York", 10, L2, bold, nohscroll, whiteBox;
	fontset = text, "New York", 12, L2, nohscroll;
	fontset = smalltext, "New York", 10, L2, nohscroll;
	fontset = input, "Courier", 12, L2, bold, nowordwrap;
	fontset = output, "Courier", 12, L2, nowordwrap;
	fontset = message, "Courier", 12, L2, R65535, nowordwrap;
	fontset = print, "Courier", 12, L2, nowordwrap;
	fontset = info, "Courier", 12, L2, nowordwrap;
	fontset = postscript, "Courier", 12, L2, nowordwrap;
	fontset = name, "Geneva", 10, L2, italic, B65535, nowordwrap, nohscroll;
	fontset = header, "Times", 10, L2;
	fontset = footer, "Times", 12, L2, center;
	fontset = help, "Geneva", 10, L2, nohscroll;
	fontset = clipboard, "New York", 12, L2;
	fontset = completions, "New York", 12, L2, nowordwrap;
	fontset = network, "Courier", 10, L2, nowordwrap;
	fontset = graphlabel, "Courier", 12, L2, nowordwrap;
	fontset = special1, "New York", 12, L2, nowordwrap;
	fontset = special2, "New York", 12, L2, center, nowordwrap;
	fontset = special3, "New York", 12, L2, right, nowordwrap;
	fontset = special4, "New York", 12, L2, nowordwrap;
	fontset = special5, "New York", 12, L2, nowordwrap;]
:[font = input; initialization; ]
*)
(* The DiracDelta and HeavisideTheta functions.
   Should be working but not guaranteed!
   PD software 3.4.1990 by
             Pekka Janhunen
             Finnish Meteorological Institute
             Geophysics Dept.
             pjanhune@finsun.csc.fi   *)

BeginPackage["DeltaTheta`"]

HeavisideTheta::usage = "HeavisideTheta[x] = 1 for x > 0 and 0 for x < 0"

DiracDelta::usage = "DiracDelta[x] = Infinity for x = 0 and 0 otherwise"
  
FourierCoefficient::usage = "FourierCoefficient[f,n,x,a,b] gives the nth
  Fourier coefficient of f when f is expanded to complex Fourier series
  in the interval [a,b]. a and b are optional, defaults are 0 and 2 Pi."
  
Begin["`Private`"]

HeavisideTheta[x_?Positive] = 1
HeavisideTheta[x_?Negative] = 0
HeavisideTheta[0] = 1/2
HeavisideTheta /:
  Limit[HeavisideTheta[f_], x_ -> a_]:= HeavisideTheta[Limit[f, x->a]]
HeavisideTheta /:
  HeavisideTheta[f_]^a_?((NumberQ[#] && Positive[#])&):= HeavisideTheta[f]

DiracDelta[x_?Positive] = 0
DiracDelta[x_?Negative] = 0
DiracDelta[0] = Infinity
DiracDelta /:
  f_[DiracDelta[x_]]:= f[0] /; x != 0

Derivative[n_][DiracDelta][x_]:= 0 /; x != 0

DiracDelta /: Limit[DiracDelta[f_], x_ -> a_]:= 0 /; Limit[f,x->a] != 0
DiracDelta /: Limit[DiracDelta[f_], x_ -> 0]:= Infinity /; Limit[f,x->a] == 0
DiracDelta /: Limit[DiracDelta[f_], x_ -> a_]:= 0 /;  Limit[f,x->a] == Infinity
DiracDelta /: Limit[DiracDelta[f_], x_ -> a_]:= 0 /;  Limit[f,x->a] == - Infinity

Unprotect[Limit]
Limit[Derivative[n_?((IntegerQ[#]&&Positive[#])&)][DiracDelta][f_], x_ -> a_]:= 0 /;
   Limit[f,x->a] != 0
Limit[Derivative[n_?((IntegerQ[#]&&Positive[#])&)][DiracDelta][f_], x_ -> 0]:= Indeterminate /; 
   Limit[f,x->a] == 0
Limit[Derivative[n_?((IntegerQ[#]&&Positive[#])&)][DiracDelta][f_], x_ -> a_]:= 0 /;
   Limit[f,x->a] == Infinity
Limit[Derivative[n_?((IntegerQ[#]&&Positive[#])&)][DiracDelta][f_], x_ -> a_]:= 0 /;
   Limit[f,x->a] == - Infinity
Protect[Limit]

Unprotect[Integrate]

(* General rule for definite integrals: *)

Integrate[f_, {x_, a_, b_}]:= 
  Block[{F = Integrate[f,x]},
        Limit[F, x->b] - Limit[F, x->a] ] /; !TrueQ[$PiecewiseIntegrate]
(* note the following: If $PiecewiseIntegrate is True, then this rule
   is unnecessary and dangerous since then indefinite integration is
   automatically called when definite integration cannot be performed
   by other means. If $PiecewiseIntegrate is False or undefined, then
   Mathematica can only do definite integrals according to this rule. *)


(* General rules for indefinite integrals of theta and delta functions: *)

Integrate[f_. HeavisideTheta[u_] + g_ HeavisideTheta[v_], x_]:=
  Integrate[f HeavisideTheta[u],x] + Integrate[g HeavisideTheta[v],x]
  
Integrate[f_.DiracDelta[u_] + g_ DiracDelta[v_], x_]:=
  Integrate[f DiracDelta[u],x] + Integrate[g DiracDelta[v],x]

Integrate[f_.HeavisideTheta[u_] + g_.DiracDelta[v_], x_]:=
  Integrate[f HeavisideTheta[u],x] + Integrate[g DiracDelta[v],x]

(* Integrals of theta functions with linear arguments: *)
Integrate[f_. HeavisideTheta[a_. x_ + b_.], x_?AtomQ]:=
  Block[{F = Integrate[f,x]},
        (F - Limit[F, x->-b/a]) HeavisideTheta[a x + b] ]

NoneZeroQ[x_List]:= Apply[And,Map[!TrueQ[#==0]&,x]]

(* Integrals of delta functions with arbitrary arguments: *)
Integrate[f_. DiracDelta[phi_], x_?AtomQ]:=
  Apply[Plus, Map[(HeavisideTheta[x-(x/.#)](f/.#)/Abs[DeltaTheta`Dphi/.#])&,
                  Solve[phi == 0, x]]]  /;
    Head[DeltaTheta`xSolved = Solve[phi == 0, x]] == List && 
    NoneZeroQ[(DeltaTheta`Dphi = D[phi,x])/. DeltaTheta`xSolved]

(* Integrals of arbitrary derivatives of delta function with linear argument: *)
Integrate[f_. Derivative[n_?((IntegerQ[#]&&Positive[#])&)][DiracDelta][a_. x_ + b_.], x_?AtomQ]:=
  Limit[f, x-> -b/a] Derivative[n-1][DiracDelta][a x + b]/a -
  (1/a) Integrate[D[f,x] Derivative[n-1][DiracDelta][a x + b], x]

Protect[Integrate]

Unprotect[Derivative]; HeavisideTheta'[x_]:= DiracDelta[x]; Protect[Derivative]

FourierCoefficient[f_, n_, x_, aa_:0, bb_:(2 Pi)]:=
  Integrate[Expand[f Exp[- 2 Pi I n x/(bb - aa)]], {x, aa, bb}]/(bb-aa)

End[]

EndPackage[]

Null


(*
^*)