(*ciste kvuli opakovanemu spousteni notebooku*)
$HistoryLength=3;
(*zakladni konstanty, asi muzou byt fixne*)
g=9.81;
km=1000.;
mm=0.001;
kW=1000.;
Ah=3600.;
(*konstanty ve vypoctu, budou je zadavat*)
mRed=15000;
ShydrOdpor=4;
Dchar=2;
cx=0.65;
ramenoValivehoOdporu=2.5mm;
rKola=300mm;
teplota=20;
aTry=0.15;
aplus=aTry;
aminus=-aTry;
(*hustota vzduchu v zavislosti na teplote*)
ro[T_]:=(341./T);

(*zadavani "jak kde chteji jet rychle", cil je mit funkci pozadovanarychlost=vz=f(x=poloha), budou zadavat v Excelu asi dva sloupecky, polohu znacime x*)
vmax=(35.5/3.6);
datRychlost={{0,0.2},{600,0.3},{1500,0.5},{3700,0.8},{3800,0.2},{8500,0.9},{11000,1},{13000,0.2},{14500,0.3},{25000,0.3}}/.{a_,b_}:>{a,b*vmax};
vz=Interpolation[datRychlost,InterpolationOrder->0];(*vz je tady funkce polohy, tedy vz=vz(x), zobrazeni o kousek niz, je po usecich konstantni, to dela InterpolationOrder->0*)

(*zdavani nadmorske vysky, vysku znacim h, zde ve tvaru {x=poloha,h=vyska}, polohu znacime x, vysku h, budou zadavat v Excelu asi dva sloupecky*)
datVyska={{0,200},{500,250},{1000,220},{3000,280},{3200,270},{5000,400},{8000,420},{9000,270},{11000,400},{12500,380},{15000,400},{25000,400}};
h=Interpolation[datVyska,InterpolationOrder->1];
(*h je tady funkce polohy, tedy h=h(x), zobrazeni o kousek niz, je po usecich primkova, to dela InterpolationOrder->1*)

(*kvuli numerickym metodam, ktere pocitaji "dopredu", beru ze posledni bod  (posledni polohu, x) je dodany jen formalne, xmax se bere z predposlednich zadanych bodu... neni to moc chytre, udelej to jakkoli nebo matlabi ODE nebude chodit dal a muze se vzit zadany konec, vymyslis to lip*)
xmax=Union@Flatten[Take[#,{-2}]&/@{datRychlost,datVyska}]//Last;

(*je videt, jaky jsem udelal trik: pozadovanou rychlost mam jako interpolaci  nulteho radu, tedy jako konstantni po usecich, kdezto nadmorskou vysku jako interpolaci prvniho radu, takze to dela kopce, je celkem jedno, jak to v Simulinku udelas, ale asi to umi interpolace uplne stejne:-)*)
plzad=Plot[vz[x],{x,0,xmax},AxesOrigin->{0,0},AxesLabel->{"x [m]","v [m.(s^(-1))]"},PlotStyle->{Blue,Thickness[0.006]}];
plh=Plot[h[x],{x,0,xmax},AxesOrigin->{0,0},AxesLabel->{"x [m]","h [m]"},PlotStyle->{Blue,Thickness[0.006]}];

(*pomocna funkce pro vypocet kinematiky, jde jen o udelani "hladkeho If:-)", v podstate jde o to, aby je-li rychlost (v) mensi, nez rychlost pozadovana (vzadana), je zrychleni az =aplus, v opacnem az =aminus, tedy vlastne regulator:-) klidne to udelej ifem, nebo jak se to bude libit ODE*)
ClearAll[vyr];
vyr[x_]:=c+b*x;
resbc=Solve[{vyr[1]==aplus,vyr[-1]==aminus}][[1]];
ka=10;
az[v_,vzadana_]:=c+b*Tanh[ka*(vzadana-v)]/.resbc;

(*jenom ukazka jak to vypada... celej spas delame, protoze skokova zmena rychlosti je nefyzikalni a tedy se musi nutne pridat nejake fyzikalne mozne pravidlo, jak menit rychlost podle nefyzikalniho pozadavku. Zde je to reseno zavedenim max kladneho a max zaporneho zrychleni... ale to neni uplne mimo, koukali jsme se Zdenkem Mullerem a Honzou Spetlikem na busy ve Svedsku a nakonec souhlasili:-), ukazano pro vzadana=50(ms^(-1))*)
Plot[az[v,50],{v,40,60},AxesLabel->{"v [m.(s^(-1))]","a [m.(s^(-2))]"},PlotStyle->{Blue,Thickness[0.006]}];

(*podstata triku:-) prijali jsme ze zrychleni je funkci aktualni a pozadovane rychlosti a napsali jsme kinematicke rovnice*)
rceRychlost={v'[t]==az[v[t],vz[x[t]]],v[0]==0,x'[t]==v[t],x[0]==0};

(*potrebujeme konec casoveho intervalu, vypocitame ho "jako kdybychom jeli skutecne pozadovanou rychlosti", o kod nejde, je to vlastne po usecich vypocitane "t=x/v" a secteno*)
tmax=Plus@@Drop[Partition[datRychlost,2,1]/.{{x1_,v1_},{x2_,v2_}}:>((x2-x1)/v1),-1];

(*jen vyreseni kinematickych difosek,MaxStepSize->0.0001tmax je jen aby num. metoda nelezla z oboru iksu, kde je definovana pozadovana rychlost vz*)
res=NDSolve[rceRychlost,{x[t],x'[t],v[t],v'[t]},{t,0,tmax},SolveDelayed->True,MaxStepSize->0.0001tmax,MaxSteps->(10^6)][[1]];

(*vypocteme, jak dlouho jsme jeli, kdyz jsme jeli, jak jsme jeli, abychom dojeli do xmax,tedy s danou zavislosti zrychleni na rychlosti pozadovane a rychlosti skutecne z kinematickych difosek*)
tm=t/.FindRoot[xmax==x[t]/.res,{t,tmax-100}];

(*ukazka jaka byla pozadovana a skutecna rychlost v terminech, jak je zadavali a asi budou takovy grafy chtit, asi i dost jinych:-)*)
plvysa=ParametricPlot[{x[t],v[t]}/.res,{t,0,tm},AspectRatio->0.8,AxesLabel->{"x [m]","v [m.(s^(-1))]"},PlotStyle->{Black,Thickness[0.005]}];

Show[plzad,plvysa,PlotRange->{{0,xmax},Automatic}];

(*vypocet vykonu (po rade) vykonu odporu vzduchu, potrebneho na jizdu do kopce, vykonu pro zrychleni, vykonu treni kolecek a z toho vyjadreneho celkoveho mechanickeho vykonu potrebneho na to aby to jelo, jak vypocitala kinematika*)
PDrug[t_]:=0.5*ro[teplota]*((v[t])^2)*ShydrOdpor*cx*v[t]/.res;
Pstoupani[t_]:=mRed*g*h'[x[t]]*v[t]/.res;
Psetrvacny[t_]:=mRed*v'[t]*v[t]/.res;
Ptreci[t_]:=mRed*g*(ramenoValivehoOdporu/rKola)*(1/(Sqrt[1+((h'[x[t]])^2)]))*v[t]/.res;
Plot[Evaluate[0.001*{PDrug[t],Pstoupani[t],Psetrvacny[t],Ptreci[t]}],{t,0,tm},PlotRange->All,AxesLabel->{"t [s]","P [kW]"}];
Psum[t_]=Plus@@{PDrug[t],Pstoupani[t],Psetrvacny[t],Ptreci[t]};
Plot[Psum[t],{t,0,tm},PlotRange->All,AxesLabel->{"t [s]","P [W"}];

plvys1=ParametricPlot[{x[t],h[x[t]]}/.res,{t,0,tm},AspectRatio->0.8,AxesLabel->{"x [m]","h [m]"},PlotStyle->{Black,Thickness[0.005]},PlotRange->All];

(* Jelikoz mame pocitat i skutecny prubeh drahy, kdyz nemame dost vykonu abychom splnili idealni prubeh a chceme vyuzit co jsme uz naprogramovali, udelame si funkci psum[x,v], coz je mechanicky vykon potrebny, kdyz v nasem bude zadanym terrainem jedeme v miste x rychlosti v se zrychlenim a, nezavisle na case t :-) 
Jednotlive vykony znaceny stejne, ale s malym pocatecnim pismenem.
Dynamiku pri reseni budeme proste psat pomoci vykonu a ne sil, jen budeme muset zacinat s nenulovou rychlosti, nebo by sily byly nekonecne. Ale to nejak obelstime:-)
psetrvacny[x,v] znamena pozadavek na zrychleni, nikoli skutecny setrvacny vykon.
*)

pDrug[v_]:=0.5*ro[teplota]*(v^2)*ShydrOdpor*cx*v;
pstoupani[x_,v_]:=mRed*g*h'[x]*v;
(*pouziti az[v,vz[x]] je dulezite pro dalsi vsechna reseni: je to vlastne P regulator zajistujici, ze dojedeme na zadanou rychlost (bude-li dost casu)*)
psetrvacny[x_,v_]:=mRed*az[v,vz[x]]*v;
ptreci[x_,v_]:=mRed*g*(ramenoValivehoOdporu/rKola)*(1/(Sqrt[1+((h'[x])^2)]))*v;
psum[x_,v_]:=pDrug[v]+pstoupani[x,v]+psetrvacny[x,v]+ptreci[x,v];


(*jen kontroly, ze to funguje a jak to funguje, v Simulinku to nebude*)
Plot[pstoupani[x[t],v[t]]/.res,{t,0,tm},PlotRange->All,AxesLabel->{"x [m]","P [W]"}];
Plot[ptreci[x[t],v[t]]/.res,{t,0,tm},PlotRange->All,AxesLabel->{"x [m]","P [W]"}];
ParametricPlot[{x[t],psetrvacny[x[t],v[t]]}/.res,{t,0,tm},PlotRange->All,AxesLabel->{"x [m]","P [W]"}];
plvys=ParametricPlot[{x[t],Psum[t]}/.res,{t,0,tm},AspectRatio->0.8,AxesLabel->{"x [m]","P [W]"},PlotStyle->{Black,Thickness[0.005]},PlotRange->All];
Show[ParametricPlot[{x[t],psum[x[t],v[t]]}/.res,{t,0,tm},PlotRange->All,AxesLabel->{"x [m]","P [W]"},AspectRatio->0.5],plvys];

(*Tu dynamiku pohybu si radeji nacvicim pomoci NDSolve, kde to vic umim, abych to pak do stavoveho modelu napsal a vedel, ze delam chybu tam v programovani a ne ve fyzice:-)
v00 je minimalni pocatecni rychlost abychom nedelili nulou pri vypoctu v'[t]*)

(*Ukazka toho, jak zavedena kinematika drzi a nedrzi rychlost...*)
(*Ve skutecnosti pomerne delikatni a nesamozrejmy trik: tim, ze je pomoci Tanh zavedeno zrychleni a setrvacny vykon jako funkce v konecnem dusledku realizujici regulator na zadanou hodnotu, muzeme simulovat fakt, ze dame jiny vykon, napr. mensi nez vykon zajistujici dany profil trasy a rychlosti na ni, presto funkce psum[x[t],v[t]] obcas alespon aplikovana zaruci snahu o zmineny profil. Snad to bude fungovat*)

(*Vykon ktery dodavam, zmenseny o pZtraty zpusobuje zrychleni*)
pZtraty[x_,v_]:=ptreci[x,v]+pstoupani[x,v]+pDrug[v];
zrychleni[pzrychlujici_,x_,v_]:=(1/(mRed*v))*(pzrychlujici-pZtraty[x,v]);

(*Tedy bude ve vysledku pzychlujici bud spravny (pokud na to mam) a nebo jiny vykon; konstrukce psum[x,t] ale zajisti snahu o dosazeni zadane rychlosti v danem miste*)
kPStiham=0.93;
rce={v'[t]==zrychleni[kPStiham*psum[x[t],v[t]],x[t],v[t]],x'[t]==v[t]};
x0=x[t]/.res/.t->0;
v00=0.1;
v0=v00+v[t]/.res/.t->0;
resDyn=NDSolve[Union[rce,{x[0]==x0,v[0]==v0}],{x[t],v[t],v'[t]},{t,0,tm},MaxSteps->(10^6),SolveDelayed->True,StartingStepSize->(10^(-9)),MaxStepSize->(10^(-4))tm];

plvysb=ParametricPlot[{x[t],v[t]}/.resDyn[[1]],{t,0,tm},AspectRatio->0.8,AxesLabel->{"x [m]","v [m.(s^(-1))]"},PlotStyle->{Magenta,Thickness[0.005]}];
Show[plvysb,plzad,plvysa,PlotRange->All];

(*budou zadavat, ale ne moc casto, jen kdyby se neco menilo na buse*)
c1Aku=42Ah;
uNom1Aku=12.4;
cSupCap=2*17.8;
RiSupCap=0.065;
UmaxSupCap=680.;
nAku=24;(*mineno serione v jedne z nParalenichAku vetvi*)
PFCmax=100.kW;
nParalenichAku=2;
PdobMax=15kW;
Post=500.;(*budou mozna chtit zadavat jako funkci, zatim neekli*)

(* celkovy maximalni naboj ve vsech aku*)

Qmax=nParalenichAku*c1Aku;

(*definice rizeni FC. Asi to nejlip povi posledni Plot*)
UFCStart=0.9*UmaxSupCap;
PFCmin=0.2PFCmax;
UFCmax=0.6UmaxSupCap;
datUcPfc={{-500,PFCmax},{UFCmax,PFCmax},{UFCStart,PFCmin},{UFCStart+1,0},{1000,0}};
pFC=Interpolation[datUcPfc,InterpolationOrder->1];
Plot[{pFC[uc],PFCmax,PFCmin},{uc,0,UmaxSupCap},PlotRange->All,PlotStyle->{{Thick,Black},{Blue,Opacity[0.3]},Red},AxesLabel->{"uc (V)","PFC (W)"},PlotLabel->"Strategie FC",Ticks->{{{UFCmax,"UFCmax"},{UFCStart,"UFCStart"}},{{PFCmax,"PFCmax"},{PFCmin,"PFCmin"}}},GridLines->{{UFCmax},None}];

(*pro vypocty el. obvodu potrebujeme znat zavislost vnitrniho odporu aku*, funkce Ri1Aku[q_] a U1Aku[q_] jsou tady jen vymyslene, abych mohl neco delat nez mi to posle v yrobce , nebo nez to nejak dostanu z grafu, muzou to byt i treba interpolacni funkce nebo podobny hnusy*)
Ri1Aku[q_]:=0.015*Exp[-2.5*((q-Qmax)/Qmax)];
RiAku[q_]:=(nAku/nParalenichAku)*Ri1Aku[q];
Plot[RiAku[q],{q,0,Qmax},PlotRange->All];
(*U1Aku[q_]:=12.4*(Tanh[4.5*(q/Qmax)]);*)
(*Musim jeste opravit, proklad jsem asi udelal pro trosku jine Qmax, ale mozna je to v poho, zatim to tam dej*)
U1Aku[q_]:=11.85+0.19* (E^(0.00014* (-144000+(q/nParalenichAku))))+((1.483* (q/nParalenichAku))/(10873+(q/nParalenichAku)));
UAku[q_]:=nAku*U1Aku[q];
Plot[UAku[q],{q,0,Qmax},PlotRange->All];

(*napeti prechodu diody, zada se jednou, dioda se v dalsim reprezentuje bud zdrojem napeti, nebo vypnuty vypinac*)
udd=1;

(*uvazovane zapojeni podle obr. 1, tak jak je to zapojeno. funkce dqduc si bere list {naboj aku, napeti na kondu, dobijeci vykon aku,vykon ze sbernice}. funkce vraci {dq,duc}, coz nejsou diferencilay, ale derivace naboje a kondu podle casu*)

(*jelikoz diodu uvazujeme jako zdroj napeti a vypinac, mame dve nahradni schemata, pro stav "rozpojeno"  obr. 2 a pro stav "spojeno" obr. 3*)

(*je to "trochu spatne, vykon pocitame z napeti bez ubytku na risc, uvazujeme ho ve stavu se spojenou diodou maly"*)
(*tahle cast je ciste o elektrice, takze se jinou dynamikou nemeni*)

ClearAll[dqduc];
dqduc[{q_,uc_,pdob_,psbernice_}]:=Module[{id,deruc,derq,isbernice,ic,riq,uAq,uq,idob,idrozpojeno,derqRozp,derucrozp,idspojeno,iAspojeno,derqSpoj,derucSpoj,spojeno,dq,duc,risc=RiSupCap,icspoj},
(*ze zadaneho vykonu ze sbernice vypoctu proud isbernice*)isbernice=((uc-(Sqrt[-4 psbernice risc+(uc^2)]))/(2 risc));

(*napoctu si parametry nahradniho schematu akumulatoru*)uq=UAku[q];riq=RiAku[q];

(*vypocteme dobijeci proud idrozpojeno z podminky:pdob==idrozpojeno×(uq+riq×idrozpojeno)*)
idrozpojeno=((-uq+(Sqrt[4 pdob riq+(uq^2)]))/(2 riq));
(*vypocteme derivaci naboje v akumulatoru a derivaci napeti na kondenzatoru*)
derqRozp=idrozpojeno;derucrozp=((-1)/cSupCap)*isbernice;
(*pro stav spojeno jsou rovnice kupodivu jednodussi:-), dobijeci proud spocteme z podminky pdob==(udd+uc)*idspojeno*)
idspojeno=(pdob/(uc+udd)); 
(*iAspojeno napeti na riq*)
iAspojeno=((-idspojeno risc+isbernice risc-uc-udd+uq)/(riq+risc));
(*derivace naboje, minus protoze proud tece "ven"*)
derqSpoj=-iAspojeno;
(*do kondu tece proud iAspojeno+idspojeno-isbernice a odtud derivace napeti kondu ve stavu "spojeno"*)
icspoj=(iAspojeno+idspojeno-isbernice);
derucSpoj=(1/cSupCap)*icspoj;
(*tady je uz jen reseno vetveni, podminka spojeni je na diode vetsi nez udd*)
spojeno=uq>uc+udd;dq=If[spojeno,derqSpoj,derqRozp];duc=If[spojeno,derucSpoj,derucrozp];{dq,duc}];
(*zkouska jak to funguje*)
dqduc[{0.6Qmax,290,15000.,20000}]

(*definice ucinnosti palivoveho clanku, bude jina, mozna i jako interpolace z tabulaky*)
(*Tady se taky nic nemeni*)
ηFC[pfc_]:=0.5;
Plot[ηFC[pfc],{pfc,0,PFCmax},AxesLabel->{"PFC (W)","ηPFC (-)"},PlotStyle->Thick,AxesOrigin->{0,0}];
(*definice icinnosti menice pro dobijeni baterii, bude jina, mozna i jako interpolace z tabulaky*)
ηdob[pdob_]:=0.95;
Plot[ηdob[pfc],{pfc,0,PdobMax},AxesLabel->{"Pdob (W)","ηdob (-)"},PlotStyle->Thick,AxesOrigin->{0,0}];
(*definice ucinnosti soustavy motor+menic, Pmech max zavedeno jen aby se mi to lip psalo, If rika, ze muze byt jina pro rekuperaci a jina pro smer vykonu ven, bude jina, mozna i jako interpolace z tabulaky*)
ηM[pmech_]:=0.9;
(*strategie dobijeni akumulatoru, tady se dobiji stale, pokud neni aku plna, mozna bude trochu jina, ale nemela by byt fci vice promennych*)
Pdob[q_,uc_]:=Module[{udobMax},udobMax=UAku[Qmax];If[And[q<Qmax,uc>udobMax],PdobMax,0]];

(*definice vykonu a zapojeni podle obr. 4*)
(*funkce derstav si bere dvojici {t=cas,stavovy vektor} a vraci derivaci stavoveho vektoru podle casu, q naboj v aku, uc napeti kondenzatoru, dalsi promenne jsou energie: Wsb energie do sbernice (spojena s kondem, vykon psbernice)),Wdob energie na dobijeni aku (vystup z menice do aku), Wfc energie z palivoveho clanku, Wzmar energie marena v odporu, kdyz bych mohl rekuperovat, ale baterka a kond vic nepoberou, Wpoh je energie do soustavy motor+menic,Wdobat je energie do nabijeciho menice, Wost energie spotrbovana na ostatni spotrebu (vse krome pohonu), Wmech je energie mechanicka pohonu, WdoFC je primarni energie vstoupivsi do palivoveho clanku *)
(*Tahle funkce sse samosebou zmeni znacne... v parametrech pribude poloha x a rychlost v, ty se budou pocitat z dynamiky*)

(*nejprve ovsem zkusime nechat soucasny stav, kdy vykon zatim staci a jen si vypocitame kinematiku*)

(*Pozor na Pmech, mel by byt fci[x,v] a nikoli fci[t] jako Psum[t], ale kvuli ladeni zatim je tak*)
(*zadani bezrozmerneho stavu napeti na kondu a stavu nabiti aku, budou zadavat, novim, jestli budou chtit takhle, nebo zadavat uc0 primo*)
socuc0=1;
soc0=1;

(*vypocteni startovaciho stavu, energie na zacatku brany nuklove, kvuli testovani derstav presunuto*)
uc0=socuc0*UmaxSupCap;
q0=Qmax*soc0;
stavStart={x0,v0,q0,uc0,0,0,0,0,0,0,0,0,0};

(*Pridavame vykonova omezeni*)
PrikonPohonMax=750*(10^3);
(*naboj v baterkach kdy jeste se snazim pozadovany vykon vytahnout z baterek*)
Qmin=0.1Qmax;
(*u maximalniho vykonu mozneho bych jiz nemel cerpat, v dalsim kdoku dostanu komplexni cisla*)
KmaxVykon=0.8;
(*PrikonPohonMax budou zadavat*)
ClearAll[derstav];
derstav[{t_,{x_,v_,q_,uc_,Wsb_,Wdob_,Wfc_,Wzmar_,Wpoh_,Wdobat_,Wost_,Wmech_,WdoFC_}}]:=Module[{Pmech,pdob,ppoh,pdobaterii,pdoFC,podminakaMareni,pDany,pfc,psbernice,pmarim,podminkaSamaSbernice,aZrychleni,pMechSkutecny,taut,uq,riq,risc,PelMax,PmechPozadovany,PelPohonPozadovany,PmaxMuzu},
(*abychom zjistili, jestli system "stiha" dodavat potrebny vykon na jizdu podle harmonogramu, napocteme maximalni dostupny elektricky vykon z vykonoveho prizpusobeni*)
uq=UAku[q];riq=RiAku[q];risc=RiSupCap;
PelMax=If[uc>uq-udd,KmaxVykon*((uc^2)/(4*risc)),If[q>Qmin,KmaxVykon*((((uq-udd))^2)/(4*riq)),0]]+pFC[uc]-Post;
(*spoctu pozadovany mechanicky a elektricky vykon "aby to jelo jak chci"*)
PmechPozadovany=psum[x,v];
PelPohonPozadovany=(PmechPozadovany/(ηM[PmechPozadovany]));
(*dobijeni ma svuj vlastni If ve sve definici, v tomto priblizeni ho pocitam uz bez ifu*)
pdob=Pdob[q,uc];
pdobaterii=(pdob/(ηdob[pdob]));
(*situace kdy pohon i baterky stihaji*)
PmaxMuzu=Min[PelMax,PrikonPohonMax];
If[PelPohonPozadovany≤PmaxMuzu,pMechSkutecny=PmechPozadovany;ppoh=PelPohonPozadovany];
(*situace kdy nestihame*)
If[PmaxMuzu<PelPohonPozadovany,pMechSkutecny=PmaxMuzu*ηM[PmaxMuzu];ppoh=(pMechSkutecny/(ηM[pMechSkutecny]))];
(*pDany je vysledny vykon sbernice, kdyby psbernice byla 0 a pfc 0, tedy kdyby se nemenilo napeti na kondu a palivovy clanek nejel*)
pDany=Post+ppoh+pdobaterii;
(*kdyby toto vychazelo zaporne, znamena to, ze mam prebytek vykonu, ktery nemam kam udat, podminkaMareni mi rika, jestli takovy stav mam*)
podminakaMareni=And[pDany≤0,uc≥UmaxSupCap];
(*plati-li podminakaMareni, nastavim pfc=0;psbernice=0;pmarim=pDany, pokud neplati, nastavim pmarim=0;pfc=pFC[uc];psbernice=pDany-pfc*)
If[podminakaMareni,pfc=0;psbernice=0;pmarim=pDany;,
pmarim=0;pfc=pFC[uc];psbernice=pDany-pfc];
(*z vysledneho nastaveneho pfc vypoctu pdoFC*)
pdoFC:=(pfc/(ηFC[pfc]));
aZrychleni=zrychleni[pMechSkutecny,x,v];
(*vektor, ktery funkce vraci, tedy vypoctene derivace naboje, derivace uc a po rade vypoctene prislusne vykony*)
{v,aZrychleni,Sequence@@dqduc[{q,uc,pdob,psbernice}],psbernice,pdob,pfc,pmarim,ppoh,pdobaterii,Post,pMechSkutecny,pdoFC}];
derstav[{0,{0.9xmax,10,0.5q0,0.9uc0,0,0,0,0,0,0,0,0,0}}]

(*nalezeni nejmensiho potrebneho casoveho kroku numericke metody*)
dv={-1,1}.({Min[#],Max[#]}&@(Last/@datRychlost));
Δt=0.005*(dv/(Max[Abs[{aminus,aplus}]]));

n=Round[(tm/(Δt))];
EuSt[stavTime:{t_,_List}]:=stavTime+Δt*{1,derstav[stavTime]};

(*jelikoz nedosahneme pokazde kvuli omezenim systemu vykonu odpovidajiciho idealnimu projeti trasy, pojedeme dele: z toho duvodu nekonci iteracni vypocet Eulerem po danem poctu casovych kroku, ale dosazenim cile, tedy x≥xmax*)
resEu=Timing[NestWhileList[EuSt,{0,stavStart},#[[2,1]]<xmax&]];
Print["trvalo to  ",resEu[[1]]," s"]
rv=resEu[[2]];
drahaRychlost=rv/.{t_,{x_,v_,q_,uc_,Wsb_,Wdob_,Wfc_,Wzmar_,Wpoh_,Wdobat_,Wost_,Wmech_,WdoFC_}}:>{x,v};
plxv=ListPlot[drahaRychlost,AxesOrigin->{0,0},PlotStyle->{Hue[0],PointSize[0.0035]}];

(*dal uz je to stejne*)

rv=resEu[[2]];
ucecka=rv/.{t_,{x_,v_,q_,uc_,Wsb_,Wdob_,Wfc_,Wzmar_,Wpoh_,Wdobat_,Wost_,Wmech_,WdoFC_}}:>{x,(uc/UmaxSupCap)};
pluc=ListPlot[ucecka,PlotRange->All,AxesLabel->{"t (s)","souc (-)"},PlotStyle->{Hue[0.63],PointSize[0.005]},GridLines->Automatic,AxesOrigin->{0,0},PlotRange->{{0,xmax},Automatic}];

socka=rv/.{t_,{x_,v_,q_,uc_,Wsb_,Wdob_,Wfc_,Wzmar_,Wpoh_,Wdobat_,Wost_,Wmech_,WdoFC_}}:>{t,(q/Qmax)};
plsoc=ListPlot[socka,PlotRange->All,AxesLabel->{"t","soc (-)"},PlotStyle->{Hue[0.63],PointSize[0.005]},GridLines->Automatic,PlotRange->All,AxesOrigin->{0,0}];
WdoFCcka=rv/.{t_,{x_,v_,q_,uc_,Wsb_,Wdob_,Wfc_,Wzmar_,Wpoh_,Wdobat_,Wost_,Wmech_,WdoFC_}}:>{t,WdoFC};
plsdoFC=ListPlot[WdoFCcka,PlotRange->All,AxesLabel->{"t","WdoFC (J)"},PlotStyle->{Hue[0.63],PointSize[0.005]},GridLines->Automatic,PlotRange->All,AxesOrigin->{0,0}];

drahaRychlost/.{x_,v_}:>v;
%//ListPlot[#,PlotRange->All,AxesLabel->{"i [-]","v [m.(s^(-1))]"},PlotStyle->PointSize[0.0055]]&;

Show[plxv,plzad,PlotRange->All];
pluc;
plzad;
plh;
plsoc;
Print["prum. rychlost= ",3.6*((rv[[-1,2,1]])/(rv[[-1,1]]))," km/h,  cas jizdy= ",rv[[-1,1]]/60," min, ujeta draha= ",(10^(-3))rv[[-1,2,1]]," km"]

datGr=rv/.{t_,{x_,v_,q_,uc_,Wsb_,Wdob_,Wfc_,Wzmar_,Wpoh_,Wdobat_,Wost_,Wmech_,WdoFC_}}:>{t,x,v,(q/Qmax),(uc/UmaxSupCap)};
poloha[{t_,x_,v_,qr_,ucr_}]:={x,h[x]};
polohy=poloha/@datGr;
dejMax:=Max[#/.{a_,b_}:>b]&;
hmax=dejMax@datVyska;
vmax=dejMax@datRychlost;
xvmax=datGr[[-1,2]];
imax=Length[datGr];(*zde hezke cislo*)
plv1=ListPlot[polohy,AxesLabel->{"x [m]","h [m]"},PlotStyle->{PointSize[0.004],Hue[0.63]},PlotRange->{{-200,1.1xvmax},{0,1.1hmax}}];
busBod[i_Integer]:=Graphics[{Hue[0],PointSize[0.025],Point[polohy[[i]]]}];
uk1[i_Integer]:=Show[plv1,busBod[i]];
Manipulate[uk1[i],{i,1,imax,1}]