(*Tak se na nas obvod podivame nejprve programatorsky inteligentneji a pak ho vyresime pro harmonicke buzeni metodami HUS*)
(*HUS (angl. SSS, tedy Sinusoidal Steady State) viz skripta CAO*)
(*ale to vam posleme az pozdeji, ted to mate jako ukol za body*)
(*podivejte, clovek, technik musi zadat rovnice (popise obvod) a hodnoty soucastek, pokud ale inteligentne zvoli format promennych a pokud jsou (jak je v elektrotechnice caste) pocatecni hodnoty velicin v derivaci nulove, jde tohle uz udelat automaticky*)

(*nase kopirovana klasika*)
Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];

(*protoze se muze stat, ze potrebujeme rovnice, vyrazy atd. v obecne forme (treba si je chcete vyexportovat do kodu LATEXu ), nezadefinujeme  globalni promenne jako napr. R1=1000, ale udelame si list pravidel (Rule) dosazeni*)
(*muze to byt zbytecne, ale sichr je sichr a tak je zadame tupu Real, tedy nap. R1->1000. namisto R1->1000, podstatna je ta tecka*)
(*proc Rule a ne RuleDelayed viz vyukovy text, ale az tak na tom nesejde: prave strany (to "za sipkou") jsou cisla nebo u zdrojunize konkretni vyrazy, je tedy jedno, jestli se vyhodnoti pri definici pravidla nebo az pri jeho pouziti*)
dosadHodnotySoucastek={R1->1000.,R2->1000.,R3->1000.,R4->1000.,R5->1000.,L1->50*(10^(-3)),L2->50.*(10^(-3)),C1-> 4.7*(10^(-6))};

(*obdobne pro zadane udaje o zdrojich, jasne, mohli bychom to napsat do soucastek... ale docela casto menime treba zdroje a ne soucastky a nekdy naopak... napriklad to muzeme vyuzit az budeme resit obvod metodami HUS i pomoci NDSolve, kdy ponechame dosazeni soucastek a bude jine dosazeni zdroju*)
dosadZdroje={uZ[t]->10*Tanh[10*Sin[1000.*t]],iZ[t]->0.01*Tanh[10*Sin[Pi*1000.*t]]};

(*a dame to dohromady*)
dosadVse=Union[dosadHodnotySoucastek,dosadZdroje];

(*pro zkoumave a chytre hlavicky: kdyz se drzime algoritmu ze skript, jsou nase rovnice  jinym popisem, ale ekvivalentim popisem naseho obvodu... takze podle techto rovnic dokazete nakreslit prislusne schema! Takze kdyz ztratime schema obvodu a zustae nam kod, nejsme ztraceni :-) *)

(*rovnice uzlu*)
uzel1=iL2[t]==iR5[t]+iZ[t];
uzel2=0==iUZ[t]+iL1[t]+iR1[t]+iL2[t];
uzel3=iUZ[t]+iR3[t]==iR2[t]+iC1[t];
uzel4=iR5[t]==iR3[t]+iR4[t]+0;

(*rovnice prvku obvodu*)
rezistorR1=u2[t]==R1*iR1[t];
rezistorR2=u3[t]==R2*iR2[t];
rezistorR3=u4[t]-u3[t]==R3*iR3[t];
rezistorR4=u4[t]==R4*iR4[t];
rezistorR5=u1[t]-u4[t]==R5*iR5[t];
kondenzatorC1=iC1[t]==C1*u3'[t];
civkaL1=u2[t]==L1*iL1'[t];
civkaL2=u2[t]-u1[t]==L2*iL2'[t];
zdrojU=uZ[t]==u3[t]-u2[t];

(*takze obvodove rovnice bez pocatecnich podminek*)
rovniceBezPocPodm={uzel1,uzel2,uzel3,uzel4,rezistorR1,rezistorR2,rezistorR4,civkaL1,kondenzatorC1,civkaL2,zdrojU,rezistorR3,rezistorR5};

(*a tyto rovnice ciselne*)
rceDosad=rovniceBezPocPodm/.dosadVse;

(*vyhledame si nezname*)
(*schvalne je chceme jako napr. u1 a ne u1[t]... NDSolve vrati o co si inteligentne rekneme a rekneme-li si o  vyraz u1[t], dostaneme pravidlo dosazeni u1[t]->Inter..., ale v tu chvili uz nemame pravidlo pro dosazeni za u1'[t] ... rekneme li si o u1, dostaneme nikoli vyraz, ale funkci a tedy napiseme-li u1'[t]/.NDSolve[rce_,vars_,interval_,opts___], tak se nejprve dosadi za funkci u1, ta se zderivuje a pak se zavola na parametr t *)
(*vzdycky je lehci dostat z funkce vyraz (proste ji zavolame na symbolicky parametr), nez opacne z vyrazu funkci*)
(*proto takhle*)
nezname=Union@Cases[rceDosad,a_Symbol[t]:>a,{0,∞}]

(*podobne se udelame pocatecni podminky, v elektrotechnice jsou casto nulove*)
pocPodminky=Union@Cases[rceDosad,a_'[t]:>a[0]==0,{0,∞}]

(*svazeme to do ranecku*)
rceDohromady=Union[rceDosad,pocPodminky];

(*reseni a Plot vysledku uz umime*)
tmax=0.01;
ktmax=3;
reseni=NDSolve[rceDohromady,nezname,{t,0,ktmax*tmax},StartingStepSize->(10^(-12)),MaxSteps->(10^6)][[1]];
{Plot[u4[t]/.reseni,{t,0,tmax}],Plot[u4[t]/.reseni,{t,0,ktmax*tmax}]}