(* ::Package:: *)

BeginPackage["TimeScales`"]
TimeScales::usage="Calculus with Time Scales. Coded by Helena Sofia Rodrigues (sofiarodrigues@esce.ipvc.pt)and Pedro A. F. Cruz(pedrocruz@ua.pt)"


(* ::Subsection:: *)
(*Read the time scale and see if a number is member of the Time Scale*)


Elemento[el_]:=If[NumberQ[el],{el},el[[1]]];
GetElements[ts_]:=Module[
{ge},
ge=Flatten[Map[Elemento,ts]];
ge
];
TSMemberQ[ts_,t_]:=MemberQ[ts,t]|| Length[Cases[IntervalMemberQ[ts,t],True]]>0||Length[Cases[Abs[t-#]<10^-16&/@ GetElements[ts],True]]>0;
TS::notints="Value `1` is not on time scale.";
TSPlot[ts_]:=Module[
{els},
Off[TS::notints];
ListaTS[ts]={#,If[TSMemberQ[ts,#],0,{}]}& /@ Range[Min[ts],Max[ts],0.01];
els=ListaTS[ts];
On[TS::notints];
ListPlot[els,AxesLabel->Automatic,Axes->{True,False},PlotLabel->"Time Scale"]
]


(* ::Section:: *)
(*Basic calculus*)


(* ::Subsection:: *)
(*Forward jump operator*)


ChooseSigma[ts_,pos_,n_,t_] :=Module[
{elem,solution},
(* Use this when t is an isolated point or the right point of an interval *)
If[pos==n,
solution := t,
elem := ts[[pos+1]];(*Print[ChooseSigma,elem];*)
solution := If[ NumberQ[elem],elem,elem[[1,1]]]
];
solution
];

TSsigma[ts_,t_]:=Module[
{pos,n,search,tsigma,intervalRightPoint},
n=Length[ts];
search=True;
pos=1;
While[search && pos<=n,
Which[

(* Main case: t is in some closed interval *)
IntervalMemberQ[ts[[pos]],t], (* comma *)
search=False;
intervalRightPoint = ts[[pos,1,2]];
tsigma = If[t<intervalRightPoint, t, ChooseSigma[ts,pos,n,t] ],

(* Main case: t is an isolated point *)
ts[[pos]]==t, (* comma *)
search = False;
tsigma = ChooseSigma[ts,pos,n,t] (* comma *)
];
pos= pos+1;
]; (* end while *)
If[search,
Message[TS::notints,t];
Return[{}]
];
tsigma (* function result *)
]

TSsigmaPlot[ts_]:=Module[
{els,elsts},
Off[TS::notints];
ListaSigma[ts]={#,TSsigma[ts,#]}& /@ Range[Min[ts],Max[ts],0.01];
els=GetElements[ts];
elsts=ListaSigma[ts];
On[TS::notints];
ListPlot[elsts,Ticks->{els},AxesLabel->{"time scale","\[Sigma]"},PlotLabel->"Forward jump operator"]
]


(* ::Subsection:: *)
(*Backward jump operator*)


ChooseRho[ts_,pos_,n_,t_] :=Module[
{elem,solution},
(* Use this when t is an isolated point or the left point of an interval *)
If[pos==1,
solution := t,
elem := ts[[pos-1]];(*Print[ChooseSigma,elem];*)
solution := If[ NumberQ[elem],elem,elem[[1,2]]]
];
solution
]

TSrho::notints="Value `1` is not on time scale.";
TSrho[ts_,t_]:=Module[
{pos,n,search,tsrho,intervalLeftPoint},
n=Length[ts];
search=True;
pos=1;
While[search && pos<=n,
Which[

(* Main case: t is in some closed interval *)
IntervalMemberQ[ts[[pos]],t], (* comma *)
search=False;
intervalLeftPoint = ts[[pos,1,1]];
tsrho = If[t>intervalLeftPoint, t, ChooseRho[ts,pos,n,t] ],

(* Main case: t is an isolated point *)
ts[[pos]]==t, (* comma *)
search = False;
tsrho = ChooseRho[ts,pos,n,t] (* comma *)

];
pos= pos+1;
]; (* end while *)
If[search,
Message[TS::notints,t];
Return[{}]
];
tsrho (* function result *)
]

TSrhoPlot[ts_]:=Module[
{els,elsts},
Off[TS::notints];
ListaRho[ts]={#,TSrho[ts,#]}& /@ Range[Min[ts],Max[ts],0.01];
els=GetElements[ts];
elsts=ListaRho[ts];
On[TS::notints];
ListPlot[elsts,Ticks->{els},AxesLabel->{"time scale","\[Sigma]"},PlotLabel->"Backward jump operator"]
]


(* ::Subsubsection:: *)
(*Set of points (Points classification)*)


TSPoint[ts_,t_]:=If[TSMemberQ[ts,t],
If[ TSsigma[ts,t]>t&&TSrho[ts,t]<t, Print[t, " is a isolated point"],
	If[TSsigma[ts,t]==t&& TSrho[ts,t]<t ,Print[t, " is left-scattered and right-dense"],
		If[TSsigma[ts,t]>t&&TSrho[ts,t]==t, Print[t, " is left-dense and right-scattered"], Print[t, " is a dense point"]]]],
 Message[TS::notints,t];
Return[{}]]


(* ::Subsubsection:: *)
(*Graininess function*)


TSTelevadok[ts_]:=If[Length[Cases[IntervalMemberQ[ts,Last[ts]],True]]>0,ts,Drop[ts,-1]]

TSmu::notints="Value `1` is not on \!\(\*SuperscriptBox[\"T\", \"k\"]\).";
TSmu[ts_,t_]:=If[
MemberQ[TSTelevadok[ts],t]||Length[Cases[IntervalMemberQ[TSTelevadok[ts],t],True]]>0,
TSsigma[ts,t]-t,
Message[TSmu::notints,t];Return[{}]]


(* ::Subsubsection:: *)
(*Backward graininess function*)


TSTindicek[ts_]:=If[Length[Cases[IntervalMemberQ[ts,First[ts]],True]]>0,ts,Drop[ts,1]]

TSnu::notints="Value `1` is not on \!\(\*SubscriptBox[\"T\", \"k\"]\).";
TSnu[ts_,t_]:=If[
MemberQ[TSTindicek[ts],t]||Length[Cases[IntervalMemberQ[TSTindicek[ts],t],True]]>0,
t-TSrho[ts,t],
Message[TSnu::notints,t];Return[{}]]


(* ::Section:: *)
(*Derivatives*)


(* ::Subsubsection:: *)
(*Delta and Nabla first derivatives*)


tsIMemberQ[ts_,t_]:=If[TSMemberQ[ts,t],TSsigma[ts,t]>t&&TSrho[ts,t]<t,Message[TS::notints,t];
Return[{}]];

tsSDMemberQ[ts_,t_]:=
If[TSMemberQ[ts,t],TSsigma[ts,t]==t&& TSrho[ts,t]<t,Message[TS::notints,t];
Return[{}]];

tsDSMemberQ[ts_,t_]:=
If[TSMemberQ[ts,t],TSsigma[ts,t]>t&&TSrho[ts,t]==t, Message[TS::notints,t];
Return[{}]];

tsDMemberQ[ts_,t_]:=
If[TSMemberQ[ts,t],TSsigma[ts,t]==t&&TSrho[ts,t]==t,Message[TS::notints,t];
Return[{}]];

TSDelta[ts_,f_,t_]:=
If[TSMemberQ[ts,t],
	If[tsIMemberQ[ts,t]||tsDSMemberQ[ts,t],(f[TSsigma[ts,t]]-f[t])/TSmu[ts,t],
	f'[t]],
Message[TS::notints,t];Return[{}]];

TSNabla[ts_,f_,t_]:=
If[TSMemberQ[ts,t],
	If[tsIMemberQ[ts,t]||tsSDMemberQ[ts,t],
(f[t]-f[TSrho[ts,t]])/TSnu[ts,t],
f'[t]],
Message[TS::notints,t];Return[{}]];



(* ::Subsubsection:: *)
(*Operations with Delta and Nabla Derivatives*)


TSDeltaSum[ts_,f_,g_,t_]:=TSDelta[ts,f,t]+TSDelta[ts,g,t];
TSDeltaProdConst[ts_,f_,k_,t_]:=k*TSDelta[ts,f,t];
TSDeltaProduct[ts_,f_,g_,t_]:=TSDelta[ts,f,t]*g[t]+f[TSsigma[ts,t]]*TSDelta[ts,g,t];
TSDeltaQuocient[ts_,f_,g_,t_]:=If[g[t]*g[TSsigma[ts,t]]!=0,(TSDelta[ts,f,t]*g[t]-f[t]*TSDelta[ts,g,t])/(g[t]*g[TSsigma[ts,t]]), Print["It is not possible to calculate"]];

TSNablaSum[ts_,f_,g_,t_]:=TSNabla[ts,f,t]+TSNabla[ts,g,t];
TSNablaProdConst[ts_,f_,k_,t_]:=k*TSNabla[ts,f,t];
TSNablaProduct[ts_,f_,g_,t_]:=TSNabla[ts,f,t]*g[t]+f[TSrho[ts,t]]*TSNabla[ts,g,t];
TSNablaQuocient[ts_,f_,g_,t_]:=If[g[t]*g[TSrho[ts,t]]!=0,(TSDelta[ts,f,t]*g[t]-f[t]*TSDelta[ts,g,t])/(g[t]*g[TSsigma[ts,t]]), Print["It is not possible to calculate"]];




EndPackage[]