(* ::Package:: *)

(* 
   This package was written by Stephan Trenn
   Version 1.0 (10.05.2010)
   Everybody is free to use and change this package. I appreciate any form of feedback.
   This package comes with a documention as a Mathematica Notebook file: documentation.nb
*)
BeginPackage["DpwSmooth`"];
  DpwSmooth::usage = "DpwSmooth is the container for all piecewise-smooth distributions, i.e. every piecewise-smooth distribution has the form 'DpwSmooth[...]'. The standard functions Plus, Times, Derivative can be applied on piecewise-smooth distributions.";
  SubD::usage = "SubD[f] returns the piecewise-smooth distribution f_D for a given smooth function f:R->R";
  ZeroD::usage = "ZeroD is the zero piecewise-smooth distribution, i.e. the neutral element corresponding to addition";
  OneD::usage = "OneD is the piecewise-smooth distribution corresponding to the function which is constantly one, i.e. the neutral element corresponding to multiplication";
  DiracImpulseAt::usage = "DiracImpulseAt[t] returns the Dirac Impulse at t as a piecewise-smooth distribution";
  DiracImpulse::usage = "DiracImpulse is the Dirac impulse (Dirac Delta function) at zero as a piecewies-smooth distribution";
  Restrict::usage = "Restrict[D,I] restricts the piecewise-smooth distribution to the interval I=Interval[...], Restrict accepts to the two optional arguments IntervalContainsLeft and IntervalContainsRight";
  IntervalContainsLeft::usage ="Optional Argument for Restrict: Restrict[D,I,IntervalContainsLeft->False], the default value for Restrict[...] is True";
  IntervalContainsRight::usage ="Optional Argument for Restrict: Restrict[D,I,IntervalContainsRight->True], the default value for Restrict[...] is False";
  Eval::usage = "Eval[D,t,\"Left\"],Eval[D,t,\"Right\"] or Eval[D,t,\"Imp\"] yields the corresponding evaluation of the piecewise-smooth distribution D at t";
  RegularPart::usage = "RegularPart[D] gives the regular part of the piecewise-smooth distribution D as a piecewise-smooth distribution";
  RegularPartAsPiecewise::usage = "RegularPartAsPiecewise[D] gives the functional part of the piecewise-smooth distribution D as a Piecewise[] function";
  ImpulsivePart::usage = "ImpulsivePart[D] gives the impulsive part of the piecewise-smooth distribution D as a piecewise-smooth distribution";
  ImpulsivePartAsDiracDelta::usage = "ImpulsivePartAsDiracDelta[D] gives the impulsive part of the piecewise-smooth distribution D as a sum of DiracDelta and their derivatives";
  DpwSmoothPlot::usage = "DpwSmoothPlot[D,I] plots the piecewise-smooth distribution D on the interval I, Dirac-Impulses are drawn as arrows whose width and color represent the different derivatives";

Begin["`Private`"];
Module[{t,imp,myIntervalMemberQ,myEqual,myMin,leftClosed,rightClosed,scalarMult,functionMerge,simpleFunctionMerge,functionSplit,simpleFunctionSplit,jumps,realAxis,singlePwConvert,functionEqual,diracsEqual,simpleImpulseEqual},
  realAxis = Interval[{-Infinity,Infinity}];
  SubD[f_]:=DpwSmooth[{{Function[t,Evaluate[f[t]]],realAxis}},{}];
  ZeroD = SubD[Function[t,0]];
  OneD = SubD[Function[t,1]];
  DiracImpulseAt[DirectedInfinity[_]]:=ZeroD;
  DiracImpulseAt[t1_] := DpwSmooth[{{Function[t,0],realAxis}},{{t1,{{0,1}}}}];
  DiracImpulse = DiracImpulseAt[0];
  Options[Restrict]={IntervalContainsLeft->True, IntervalContainsRight->False};

  DpwSmooth[]=ZeroD;
  DpwSmooth[f_Function]:=SubD[f];
  DpwSmooth[DiracDelta]:=DiracImpulse;
  DpwSmooth[DiracDelta[Plus[t_NumericQ,_]]]:=DiracImpulseAt[t];
  DpwSmooth[DiracDelta[Plus[_,t_NumericQ]]]:=DiracImpulseAt[t];
  DpwSmooth[DiracDelta[_]]:=DiracImpulse;
  DpwSmooth[HeavisideTheta]:=Restrict[OneD,Interval[{0,Infinity}]];

  myEqual[a_DirectedInfinity,b_] := a == b;
  myEqual[a_,b_DirectedInfinity] := a == b;
  myEqual[a_?NumericQ,b_?NumericQ] := a==b || Abs[a-b]/(1+Max[Abs[a],Abs[b]]) < 2 $MachineEpsilon;
  myEqual[i1_Interval,i2_Interval] := myEqual[Min[i1],Min[i2]] && myEqual[Max[i1],Max[i2]];
  myEqual::internalerror = "[?NumericQ,?NumericQ], [+-Infinity,_], [_,+-Infinity]  or [Interval,Interval] expected, invalid argument `1`";
  myEqual[arg___] := Message[myEqual::internalerror,{arg}];

  myMin[i_Interval] := Min[i];
  myMin::internalerror = "Interval expected, invalid argument `1`";
  myMin[arg___] := Message[myMin::internalerror,{arg}];
  
  Options[myIntervalMemberQ]={leftClosed->True, rightClosed->False};
  myIntervalMemberQ[testinterval_Interval,testt_?NumericQ,OptionsPattern[]]:=Module[{returnvalue},
    returnvalue:=If[OptionValue[leftClosed], 
      If[!OptionValue[rightClosed],
        IntervalMemberQ[testinterval,testt] && !myEqual[testt,Max[testinterval]],
        IntervalMemberQ[testinterval,testt]
	  ],
      If[!OptionValue[rightClosed],
        IntervalMemberQ[testinterval,testt] && !myEqual[testt,Min[testinterval]] && !myEqual[testt,Max[testinterval]],
        IntervalMemberQ[testinterval,testt] && !myEqual[testt,Min[testinterval]]
	  ]
    ];
	returnvalue
  ];
  myIntervalMemberQ::internalerror = "[Interval,?NumericQ,?OptionQ] expected, invalid argument `1`";
  myIntervalMemberQ[arg___] := Message[myIntervalMemberQ::internalerror,{arg}];

  simpleFunctionMerge[head_List,{fi_Function,fiInterval_Interval}]/;MatchQ[Last[head],{_Function,_Interval}]:=Module[{returnvalue,lastf,lastInterval},
    lastf = First[Last[head]];
    lastInterval = Last[Last[head]];
    returnvalue := If[Evaluate[lastf[t]]===Evaluate[fi[t]],
      Append[Most[head],{fi,Interval[{Min[lastInterval],Max[fiInterval]}]}],
      Append[head,{fi,fiInterval}]
	];
	returnvalue
  ];
  simpleFunctionMerge::internalerror= "[List,{Function,Interval}] with Last[List]={Function,Interval} expected, invalid argument `1`";
  simpleFunctionMerge[arg___] := Message[simpleFunctionMerge::internalerror,{arg}];
	
  functionMerge[f_]:=Fold[simpleFunctionMerge,{First[f]},Rest[f]];

  jumps[f_List] := Complement[Map[myMin[Last[#]]&,f],{-Infinity}];
  jumps::internalerror = "List expected, invalid argument `1`";
  jumps[arg___] := Message[jumps::internalerror,{arg}];

  simpleFunctionSplit[{fi_Function,fiInterval_Interval},t_?NumericQ]:=Module[{returnvalue},
	returnvalue:=If[myIntervalMemberQ[fiInterval,t,{leftClosed->False,rightClosed->False}],
	  {{fi,Interval[{Min[fiInterval],t}]},{fi,Interval[{t,Max[fiInterval]}]}},
	  {{fi,fiInterval}}
	];
    returnvalue
  ];
  simpleFunctionSplit::internalerror = "[{Function,Interval},NumericQ] expected, invalid argument `1`";
  simpleFunctionSplit[arg___] := Message[simpleFunctionSplit::internalerror,{arg}];

  functionSplit[f_List,t_?NumericQ]:=Flatten[Map[simpleFunctionSplit[#,t]&,f],1];

  Restrict[d_DpwSmooth,i_Interval,opts___]/;d==ZeroD := ZeroD;
  Restrict[DpwSmooth[f_List,diracs_List],Interval[i1_,rest__],opts___]:=Restrict[DpwSmooth[f,diracs],Interval[i1],opts]+Restrict[DpwSmooth[f,diracs],Interval[rest],opts];
  Restrict[DpwSmooth[f_List,diracs_List],interval_Interval,OptionsPattern[]]/;(NumericQ[(Max[interval]/.Infinity->1)]&&NumericQ[(Min[interval]/.-Infinity->-1)]) := Module[{singleFunctionRestrict,nonEmptyInterval,makeGlobal,functionRestrict,impulseRestrict,singleFun},    
    singleFunctionRestrict[{fi_Function,fiInterval_Interval}]:={fi,IntervalIntersection[fiInterval,interval]};
    singleFunctionRestrict::internalerror = "{{Function,Interval}} expected, invalid argument `1`";
    singleFunctionRestrict[arg___] := Message[singleFunctionRestrict::internalerror,{arg}];
    
    nonEmptyInterval[Interval[]] = False;
    nonEmptyInterval[Interval[{a_,a_}]] := False;
    nonEmptyInterval[Interval[{a_,Infinity}]] := a!=Infinity;
    nonEmptyInterval[Interval[{-Infinity,a_}]] := a!=-Infinity;
    nonEmptyInterval[Interval[{a_?NumericQ,b_?NumericQ}]]:=Not[myEqual[a,b]];
    nonEmptyInterval::internalerror = "Interval[{?NumericQ,?NumericQ}] expected, invalid argument `1`";
    nonEmptyInterval[arg___] := Message[nonEmptyInterval::internalerror,{arg}];
    makeGlobal[{{singleFun_Function,i_Interval}}] /; nonEmptyInterval[i] := 
      If[singleFun[t]===0,
        {{singleFun,Interval[{-Infinity,Infinity}]}},
        If[Min[i]==-Infinity,
          If[Max[i]==Infinity,
            {{singleFun,i}},
            {{singleFun,i},{Function[t,0],Interval[{Max[i],Infinity}]}}  
          ],  
          If[Max[i]==Infinity,
            {{Function[t,0],Interval[{-Infinity,Min[i]}]},{singleFun,i}},
            {{Function[t,0],Interval[{-Infinity,Min[i]}]},{singleFun,i},{Function[t,0],Interval[{Max[i],Infinity}]}}
          ]
        ]
      ];
    makeGlobal[{{leftFun_Function,leftInt_Interval},Middle___,{rightFun_,rightInt_Interval}}]:=Module[{leftPart,rightPart},
      leftPart := If[Min[leftInt]==-Infinity,
        {leftFun,leftInt},
        If[leftFun[t]===0,
          {leftFun,Interval[-Infinity,Max[leftInt]]},
          Sequence[{Function[t,0],Interval[-Infinity,Min[leftInt]]},{leftFun,leftInt}]]];
      rightPart := If[Max[rightInt]==Infinity,
        {rightFun,rightInt},
        If[rightFun[t]===0,
          {rightFun,Interval[Min[rightInt],Infinity]},
          Sequence[{rightFun,rightInt},{Function[t,0],Interval[Max[rightInt],Infinity]}]]];
      {leftPart,Middle,rightPart}
    ];
    makeGlobal::internalerror = "List of {Function,Interval} with non-empty Interval expected, invalid argument `1`";
    makeGlobal[arg___] := Message[makeGlobal::internalerror,{arg}];
    functionRestrict := makeGlobal[Select[Map[singleFunctionRestrict,f],nonEmptyInterval[Last[#]]&]];
    impulseRestrict := Select[diracs,myIntervalMemberQ[interval,First[#],{leftClosed->OptionValue[IntervalContainsLeft], rightClosed->OptionValue[IntervalContainsRight]}]&];
    DpwSmooth[functionRestrict,impulseRestrict]
  ];
  Restrict[f_Function,interval_Interval,opts___]:=Restrict[SubD[f],interval];

  RegularPart/:Plus[RegularPart[d1_],ImpulsivePart[d2_]]/;d1==d2:=d1;
  ImpulsivePart/:Plus[ImpulsivePart[d1_],RegularPart[d2_]]/;d1==d2:=d1;
    DpwSmooth/:Plus[DpwSmooth[f1_List,diracs1_List],DpwSmooth[f2_List,diracs2_List]]:=Module[{simpleFunctionAdd,impulseMerge,simpleImpulseMerge,singleImpulseMerge,simpleImpulseAdd},    	
	simpleFunctionAdd[{{fi1_Function,fi1Interval_Interval},{fi2_Function,fi2Interval_Interval}}]/;myEqual[fi1Interval,fi2Interval] :=
      {Function@@{t,Evaluate[fi1[t]+fi2[t]]},fi1Interval};
    simpleFunctionAdd::internalerror = "[{{Function1,Interval1},{Function2,Interval2}}] with Interval1==Interval2 expected, invalid argument `1`";
    simpleFunctionAdd[arg___] := Message[simpleFunctionAdd::internalerror,{arg}];
    
	simpleImpulseAdd[imps_List]:=
	  Fold[
	    If[First[Last[#1]]==First[#2],
		  Append[Most[#1],{First[#2],Last[Last[#1]]+Last[#2]}],
		  Append[#1,#2]
        ]&,
        {First[imps]},
        Rest[imps]
      ];
	
	singleImpulseMerge[imp1_List,imp2_List]:=Select[simpleImpulseAdd[Sort[Join[imp1,imp2],First[#1]<First[#2]&]],Last[#]=!=0&];

	simpleImpulseMerge[imps_List,imp_List]/;MatchQ[Last[imps],_List](*&&NumericQ[First[Last[imps]]]&&NumericQ[First[imp]]*):=
	  If[First[Last[imps]]===First[imp],
		Append[Most[imps],{First[imp],singleImpulseMerge[Last[Last[imps]],Last[imp]]}],
		Append[imps,imp]
	  ];
    simpleImpulseMerge::internalerror= "[List,{NumericQ,List}] with Last[List]={NumericQ,List} expected, invalid argument `1`";
	simpleImpulseMerge[arg___] := Message[simpleImpulseMerge::internalerror,{arg}];

    impulseMerge[{}]:={};
	impulseMerge[imps_List]:=Select[Fold[simpleImpulseMerge,{First[imps]},Rest[imps]],Last[#]!= {}&];

    DpwSmooth[
       functionMerge[Map[simpleFunctionAdd,Transpose[{Fold[functionSplit,f1,jumps[f2]],Fold[functionSplit,f2,jumps[f1]]}]]],
       impulseMerge[Sort[Join[diracs1,diracs2],First[#1]<First[#2]&]]
    ]
  ];

  scalarMult[a_,DpwSmooth[f_List,diracs_List]]:=
    DpwSmooth[
      Map[{t\[Function]Evaluate[a First[#][t]],Last[#]}&,f],
	  Map[imp\[Function]{First[imp],Map[{First[#],a Last[#]}&,Last[imp]]},diracs]
    ];

  DpwSmooth::nonCommutativeMultiplication = "Use ** or NonCommutativeMultiply[...] for multiplication of piecewise-smooth distributions";
  DpwSmooth/:Times[d1_DpwSmooth,d2_DpwSmooth]:=Module[{},
     Message[DpwSmooth::nonCommutativeMultiplication];
	 HoldForm[d1 d2]
  ];
  DpwSmooth/:Times[a_,d_DpwSmooth]:=scalarMult[a,d];
  DpwSmooth/:Times[d_DpwSmooth,a_]:=scalarMult[a,d];

  Eval[DpwSmooth[{{f_Function,realAxis}},_List],t_,"Left"]:=f[t];
  Eval[DpwSmooth[f_List,_List],t_?NumericQ,"Left"]:=Plus@@Map[
    (If[myIntervalMemberQ[Last[#],t,{leftClosed->False,rightClosed->True}],First[#][t],0])&,
    f];
  Eval[DpwSmooth[{{f_Function,realAxis}},_List],t_,"Right"]:=f[t];
  Eval[DpwSmooth[f_List,_List],t_?NumericQ,"Right"]:=Plus@@Map[
    (If[myIntervalMemberQ[Last[#],t,{leftClosed->True,rightClosed->False}],First[#][t],0])&,
    f];
  Eval[DpwSmooth[_List,diracs_List],t1_,"Imp"]:=
    DpwSmooth[{{Function[t,0],Interval[{-Infinity,Infinity}]}},Select[diracs,t1==First[#]&]];

  Derivative[1][DpwSmooth[f_List,diracs_List]]:=Module[{functionDerivative,functionImpulses,impulseDerivative,jumpTimes,jumpValues,jumpPairs},
    functionDerivative:=functionMerge[Map[{(t\[Function]Evaluate[D[First[#][t],t]]),Last[#]}&,f]];
    jumpTimes:=jumps[f];
	jumpValues:=Map[(Eval[DpwSmooth[f,{}],#,"Right"]-Eval[DpwSmooth[f,{}],#,"Left"])&,jumpTimes];
	jumpPairs:=Select[Transpose[{jumpTimes,jumpValues}],(Last[#]=!=0)&];
	functionImpulses:=Map[{First[#],{{0,Last[#]}}}&,jumpPairs];
    impulseDerivative:=Map[imp\[Function]{First[imp],Map[{First[#]+1,Last[#]}&,Last[imp]]},diracs];
    Plus[DpwSmooth[functionDerivative,functionImpulses],DpwSmooth[{{Function[t,0],realAxis}},impulseDerivative]]
  ];
  Derivative[n_Integer][d_DpwSmooth]/;(n>1):=Derivative[1][Derivative[n-1][d]];

  Derivative[1][Restrict[DpwSmooth[f_List,diracs_List],i_Interval]]:=
    Restrict[DpwSmooth[f,diracs]',i]
    + Times[Eval[DpwSmooth[f,diracs],Min[i],"Right"],DiracImpulseAt[Min[i]]]
    + Times[-Eval[DpwSmooth[f,diracs],Max[i],"Left"],DiracImpulseAt[Max[i]]];

  Derivative[n_Integer][Restrict[d_DpwSmooth,i_Interval]]/;(n>1):=Derivative[1][Derivative[n-1][Restrict[d,i]]];
 
 (* The following might conflict with other packages, however its a generic property of the derivative with respect to multiplication *)
  Derivative[1][d1_**d2_]:=d1'**d2 + d1**d2';
  Derivative[n_Integer][d1_**d2_]/;(n>1):=Derivative[n-1][d1'**d2] + Derivative[n-1][d1**d2'];

  DpwSmooth/:NonCommutativeMultiply[DpwSmooth[f1_List,diracs1_List],DpwSmooth[f2_List,diracs2_List]]:=Module[{functionMultiply,leftImpMult,rightImpMult,simpleFunctionMultiply,singleImpMult,simpleImpMult},
    
	simpleImpMult[{impOrder_,impValue_},f_,timp_,leftright_]:=Module[{ft,returnvalue},
	  ft:=Eval[DpwSmooth[f,{}],timp,leftright];
	  If[impOrder==0,
        If[NumericQ[ft]&&ft==0,
          returnvalue=ZeroD,
		  returnvalue=DpwSmooth[{{Function[t,0],realAxis}},{{timp,{{0,ft impValue}}}}]
        ],
        returnvalue = Plus[Derivative[1][simpleImpMult[{impOrder-1,impValue},f,timp,leftright]],
         simpleImpMult[{impOrder-1,-impValue},First[DpwSmooth[f,{}]'],timp,leftright]]
      ];
	  returnvalue
    ];

	singleImpMult[imp_List,f_List,leftright_]:=Fold[Plus,ZeroD,Map[(simpleImpMult[#,f,First[imp],leftright])&,Last[imp]]];
	leftImpMult[imps_List,f_List]:=Fold[Plus,ZeroD,Map[(singleImpMult[#,f,"Left"])&,imps]];
	rightImpMult[f_List,imps_List]:=Fold[Plus,ZeroD,Map[(singleImpMult[#,f,"Right"])&,imps]];
	
	simpleFunctionMultiply[{{fi1_Function,fi1Interval_Interval},{fi2_Function,fi2Interval_Interval}}]/;myEqual[fi1Interval,fi2Interval] :=
      {Function@@{t,Evaluate[fi1[t] fi2[t]]},fi1Interval};
    simpleFunctionMultiply::internalerror = "[{{Function1,Interval1},{Function2,Interval2}}] with Interval1==Interval2 expected, invalid argument `1`";
    simpleFunctionMultiply[arg___] := Message[simpleFunctionMultiply::internalerror,{arg}];

    functionMultiply[g1_,g2_]:=
      functionMerge[Map[simpleFunctionMultiply,Transpose[{Fold[functionSplit,g1,jumps[g2]],Fold[functionSplit,g2,jumps[g1]]}]]];

    Plus[DpwSmooth[functionMultiply[f1,f2],{}],leftImpMult[diracs1,f2],rightImpMult[f1,diracs2]]
  ];
  RegularPart[DpwSmooth[f_List,_]]:=DpwSmooth[f,{}];
  RegularPartAsPiecewise[DpwSmooth[f_List,_]]:=
    t\[Function]Evaluate[Piecewise[Map[{First[#][t],Min[Last[#]]<=t&&If[Max[Last[#]]==Infinity,t<=Infinity,t<Max[Last[#]]]}&,f]]];
  
  ImpulsivePart[DpwSmooth[_,diracs_List]]:=DpwSmooth[{{Function[t,0],realAxis}},diracs];
  ImpulsivePartAsDiracDelta[DpwSmooth[_,diracs_List]]:=t\[Function]Evaluate[Plus@@Map[imp\[Function]Plus@@Map[Last[#] D[DiracDelta[t-First[imp]],{t,First[#]}]&,Last[imp]],diracs]];
  
  DpwSmoothPlot[DpwSmooth[f_List,diracs_List],i_Interval,opts___]:=Module[{plotSingleImp, plotImps,myColorList,plotFunction},
    myColorList:=ColorData["Indexed","ColorList"][[34]];
    plotSingleImp[t_?NumericQ,{impOrder_Integer,impValue_?NumericQ}]:=Module[{color,thickness,arrowheads},
      color:=myColorList[[Min[impOrder+1,Length[myColorList]]]];
      thickness:=Thickness[0.005 (impOrder+1)];
      arrowheads:=Arrowheads[0.02 (impOrder+1)];
      {color,thickness,arrowheads,Arrow[{{t,0},{t,impValue}}]}
    ];
    plotImps[t_?NumericQ,imps_List]:=Map[(plotSingleImp[t,#])&,Reverse[imps]];
    plotFunction=First[Plot[RegularPartAsPiecewise[DpwSmooth[f,{}]][t],{t,Min[i],Max[i]},opts]];
    
    Graphics[{plotFunction,Map[(plotImps[First[#],Last[#]])&,diracs]},Axes->True]
  ];

  Format[DpwSmooth[{{f_Function,realAxis}},{}]]:=
    If[MatchQ[Last[f],_?NumericQ],
      If[Last[f]==0,
        StringForm["ZeroD"],
        If[Last[f]==1,
          StringForm["OneD"],
          If[Last[f]==-1,
            StringForm["-OneD"],
            StringForm["`` OneD",StandardForm[f[t]]]
          ]
        ]
      ],
      If[FreeQ[{Last[f]},First[f]],
        If[MatchQ[Last[f],_Plus],
          StringForm["(``)OneD",Last[f]],
          StringForm["`` OneD",Last[f]]],
        StringForm["SubD[t\[Function]``]",StandardForm[f[t]]/.t->StringForm["t"]]
      ]
    ];
  Format[DpwSmooth[f_List,{}]]:=Module[{restrictedFormat,addTerm,nonZeroParts},
    restrictedFormat[{fi_Function,fiInterval_Interval}]:=StringForm["Restrict[``,``]",StringForm["t\[Function]``",StandardForm[fi[t]]/.t->StringForm["t"]],StandardForm[fiInterval]];
    addTerm[oldString_StringForm,{fi_Function,fiInterval_Interval}]:=StringForm["`` + ``",oldString,restrictedFormat[{fi,fiInterval}]];
    nonZeroParts:=Select[f,Last[First[#]]=!=0&];
    Fold[addTerm,restrictedFormat[First[nonZeroParts]],Rest[nonZeroParts]]
  ];

  Format[DpwSmooth[f_List,diracs_List]]:=Module[{functionFormat,diracsFormat,maxPrimes,primesString,addSimpleTerm,simpleDiracFormat,addSingleTerm,singleDiracFormat,signString,diracsFormatFirst},
    functionFormat:=Format[DpwSmooth[f,{}]];
    maxPrimes=5;
    primesString[0]="";
    primesString[n_Integer]:=StringJoin[primesString[n-1],"'"];
    singleDiracFormat[t_,{n_Integer,value_},pos_String]:=Module[{valueString,diracString},
      valueString:=Module[{returnvalue},
        If[MatchQ[value,_Plus],
          returnvalue=StringForm["(``)",value],
          If[NumericQ[value],
            If[Abs[value]==1,
              returnvalue="",
              If[value<0,
                returnvalue=StringForm["`` ",-value],
                returnvalue=StringForm["`` ",value]
              ];
            ],
            If[MatchQ[value,Times[-1,_]],
              returnvalue=StringForm["`` ",Last[value]],
              returnvalue=StringForm["`` ",value]
            ];
          ];
        ];
        returnvalue
      ];
      signString:=If[(NumericQ[value]&&value<0) || MatchQ[value,Times[-1,_]],If[pos=="first","-","- "],If[pos=="first","","+ "]];
      diracString:=If[t===0,StringForm["DiracImpulse"],StringForm["DiracImpulseAt[``]",t]];
      If[n>maxPrimes,StringForm["````Derivative[``][``]",signString,valueString,n,diracString],StringForm["````````",signString,valueString,diracString,primesString[n]]]
    ];
    addSingleTerm[t_,oldString_StringForm,{n_Integer,value_}]:=StringForm["`` ``",oldString,singleDiracFormat[t,{n,value},"second"]];
    simpleDiracFormat[{t_,simpleDirac_List},pos_String]:=Fold[addSingleTerm[t,#1,#2]&,singleDiracFormat[t,First[simpleDirac],pos],Rest[simpleDirac]];
    addSimpleTerm[oldString_StringForm,simpleDirac_List]:=StringForm["`` ``",oldString,simpleDiracFormat[simpleDirac,"second"]];
    diracsFormat:=Fold[addSimpleTerm,simpleDiracFormat[First[diracs],"second"],Rest[diracs]];
    diracsFormatFirst:=Fold[addSimpleTerm,simpleDiracFormat[First[diracs],"first"],Rest[diracs]];
    If[Length[f]==1 && Last[First[First[f]]]===0,StringForm["``",diracsFormatFirst],StringForm["`` ``",functionFormat,diracsFormat]]
  ];
 simpleImpulseEqual[{t1_,imps1_List},{t2_,imps2_List}]:= t1==t2 && Length[imps1]==Length[imps2] && 
   And@@Map[(First[First[#]]==First[Last[#]] && Last[First[#]]==Last[Last[#]])&,Transpose[{imps1,imps2}]];
 diracsEqual[d1_List,d2_List]:=Length[d1]==Length[d2] && And@@Map[simpleImpulseEqual[First[#],Last[#]]&,Transpose[{d1,d2}]];
 functionEqual[f1_List,f2_List]:=Length[f1]==Length[f2] && And@@Map[(Last[First[#]]===Last[Last[#]] && First[First[#]][t]==First[Last[#]][t])&, Transpose[{f1,f2}]];
 DpwSmooth/:Equal[DpwSmooth[f1_List,diracs1_List],DpwSmooth[f2_List,diracs2_List]]:=functionEqual[f1,f2] && diracsEqual[diracs1,diracs2];
]; (* Private Module *)
End[];
EndPackage[];