SetOptions[EvaluationNotebook[],Background->Black,FontColor->White]

cm=0.01;
mm=0.001;
mm2=10^-6;
m=1;

Needs["SymbolicC`"]

S=758.09mm2; (*785AL1/43ST1A*)
Rv1=(2.9/2)*mm;(*poloměr vodičů v 1. a 2. vrstvě *)
Rv2=(4.12/2)*mm;(*poloměr vodičů ve 3. a 4. vrstvě*)
a=200mm;(*rozteč svazku vodorovná*)
b=200mm;(*rozteč svazku svislá*)
Dg=800m;(*hloubka fiktivního vodiče*)
Si1=6.605mm2;(*průřez vodičů v 1. a 2. vrstvě*)
Si2=13.332mm2;(*průřez vodičů v 1. a 2. vrstvě*)
Rho=0.0285*10^-6;(*ρ hliníku*)
Sal=758.09*mm2;(*průřez Al vodičů *)
Ksi=0.81;
f=50;
ω=2.Pi*f;

TypSvazku=0;(*Symetrický svazek = 1, Nesymetrický svazek = 0*)
PocetLanVeSvazku=4;

(*1. vrstva A*)

x1[i_]=5.66mm*Cos[((Pi*i)/6)+(Pi/12)];
y1[i_]=5.66mm*Sin[((Pi*i)/6)+(Pi/12)];(*pozice jednotlivých drátů v laně*)

Xvrstva1=Table[x1[i],{i,0,11}];
Yvrstva1=Table[y1[i],{i,0,11}];

Vrstva1=Transpose@{Xvrstva1,Yvrstva1};(*pozice jednotlivých drátů v laně - seznam souřadnic - 1. vrstva*)
ListPlot[Vrstva1,AspectRatio->1]

(*2. vrstva A*)
x2[i_]=8.56mm*Cos[((Pi*i)/9)];
y2[i_]=8.56mm*Sin[((Pi*i)/9)];

Xvrstva2=Table[x2[i],{i,0,17}];
Yvrstva2=Table[y2[i],{i,0,17}];

Vrstva2=Transpose@{Xvrstva2,Yvrstva2};(*pozice jednotlivých drátů v laně - seznam souřadnic - 2. vrstva*)
(*ListPlot[{Vrstva1,Vrstva2},AspectRatio->1]*)

(*3. vrstva A*)
x3[i_]=12.07mm*Cos[((Pi*i)/9)+(Pi/18)];
y3[i_]=12.07mm*Sin[((Pi*i)/9)+(Pi/18)];

Xvrstva3=Table[x3[i],{i,0,17}];
Yvrstva3=Table[y3[i],{i,0,17}];

Vrstva3=Transpose@{Xvrstva3,Yvrstva3};(*pozice jednotlivých drátů v laně - seznam souřadnic - 3. vrstva*)
(*ListPlot[{Vrstva1,Vrstva2,Vrstva3},AspectRatio->1]*)

(*4. vrstva A*)
x4[i_]=16.19mm*Cos[((Pi*i)/12)];
y4[i_]=16.19mm*Sin[((Pi*i)/12)];

Xvrstva4=Table[x4[i],{i,0,23}];
Yvrstva4=Table[y4[i],{i,0,23}];

Vrstva4=Transpose@{Xvrstva4,Yvrstva4};(*pozice jednotlivých drátů v laně - seznam souřadnic - 4. vrstva*)
ListPlot[{Vrstva1,Vrstva2,Vrstva3,Vrstva4},AspectRatio->1]
LanoAlFe=lanoA=Join[Vrstva1,Vrstva2,Vrstva3,Vrstva4];

(*1-svazek symetrický*)
svazek1S=LanoAlFe;

(*2-svazek symetrický*)
VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{-0.5*a,0}
]]]][[2,1]];
lanoA1=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{0.5*a,0}
]]]][[2,1]];
lanoA2=VypLana/.Null->Sequence[];

svazek2S=Join[lanoA1,lanoA2];


(*3-svazek symetrický*)
VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{0,0.5*b}
]]]][[2,1]];
lanoA1=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{-0.5*a,-0.5*b}
]]]][[2,1]];
lanoA2=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{0.5*a,-0.5*b}
]]]][[2,1]];
lanoA3=VypLana/.Null->Sequence[];

svazek3S=Join[lanoA1,lanoA2,lanoA3];


(*4-svazek symetrický*)
VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{-0.5*a,0.5*b}
]]]][[2,1]];
lanoA1=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{0.5*a,0.5*b}
]]]][[2,1]];
lanoA2=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{-0.5*a,-0.5*b}
]]]][[2,1]];
lanoA3=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{0.5*a,-0.5*b}
]]]][[2,1]];
lanoA4=VypLana/.Null->Sequence[];

svazek4S=Join[lanoA1,lanoA2,lanoA3,lanoA4];


(*3-svazek nesymetrický*)
VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{-a,0}
]]]][[2,1]];
lanoA1=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{0,0}
]]]][[2,1]];
lanoA2=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{a,0}
]]]][[2,1]];
lanoA3=VypLana/.Null->Sequence[];

svazek3N=Join[lanoA1,lanoA2,lanoA3];


(*4-svazek nesymetrický*)
VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{-1.5*a,0}
]]]][[2,1]];
lanoA1=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{-0.5*a,0}
]]]][[2,1]];
lanoA2=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{0.5*a,0}
]]]][[2,1]];
lanoA3=VypLana/.Null->Sequence[];

VypLana=Reap[Sow[For[i=1,i≤Length[LanoAlFe],i++,Sow[
LanoAlFe[[i]]+{1.5*a,0}
]]]][[2,1]];
lanoA4=VypLana/.Null->Sequence[];

svazek4N=Join[lanoA1,lanoA2,lanoA3,lanoA4];

svazek=If[TypSvazku==1,
If[PocetLanVeSvazku==1,svazek1S,If[PocetLanVeSvazku==2,svazek2S,If[PocetLanVeSvazku==3,svazek3S,If[PocetLanVeSvazku==4,svazek4S,"špatná hodnota PocetLanVeSvazku S"]]]],
If[PocetLanVeSvazku==3,svazek3N,If[PocetLanVeSvazku==4,svazek4N,"špatná hodnota PocetLanVeSvazku N"]]];


ListPlot[svazek,AspectRatio->1]

(*rozdíl souřadnic jednotlivých bodů*)
radkysnull=Reap[For[j=1,j≤Length[svazek],j++,Sow[For[i=1,i≤Length[svazek],i++,Sow[
svazek[[i]]-svazek[[j]]
]
]]]
][[2,1]];
radky=radkysnull/.Null->Sequence[];
Length[radky];


(*Vzdálenosti mezi jednotlivými vodici*)

PocCyklu=72*PocetLanVeSvazku;
matice=Reap[For[i=1,i≤Length[radky],i++,Sow[

(Sqrt[{radky[[i,1]]}^2+{radky[[i,2]]}^2])

]];][[2,1]];
Length[matice];
maticeVzdalenosti=Partition[matice,PocCyklu];
maticeVzdalenosti//MatrixForm;

listpromatici1=Table[Rv1,30];
listpromatici2=Table[Rv2,42];(*Příprava diagonály do matice vzdáleností matice*)
diagonalaJedenVodic=Join[listpromatici1,listpromatici2];
diagonala01=Table[diagonalaJedenVodic,PocetLanVeSvazku];
diagonala1=Flatten[diagonala01];
Length[diagonala1];


maticediagonalni=DiagonalMatrix[diagonala1];

vzdalenostniMatice=maticeVzdalenosti+maticediagonalni;
vzdalenostniMatice//MatrixForm;(*hotová matice obsahující vzdálenosti mezi všemi dráty ve svazku*)

(*Vypocet L matice*)
MatIndukcnosti=maticediagonalni;

For[j=1,j≤PocCyklu,j++,Sow[For[i=1,i≤PocCyklu,i++,

MatIndukcnosti[[i,j]]=I*0.001*ω*0.46*Log[10,(Dg/(vzdalenostniMatice[[i,j]]))]

]]];



For[i=1,i≤PocCyklu,i++,Sow[

MatIndukcnosti[[i,i]]=(Rho/(Pi*vzdalenostniMatice[[i,i]]^2))*1000+I*ω*0.001*0.46*Log[10,(Dg/(Ksi*vzdalenostniMatice[[i,i]]))]

]];
MatrixForm[MatIndukcnosti];

(*Pridani radku do L matice*)
radekdeltaU1=Table[{1},PocCyklu];
radekdeltaU=Insert[radekdeltaU1,0,PocCyklu+1];
Maticevypocet=Insert[MatIndukcnosti,radekdeltaU1,PocCyklu+1];

(*Pridani sloupce do L matice*)
sloupecMatice=Join[Table[{-1},PocCyklu],{{0}}];
Lmatice1=MapThread[Append,{Maticevypocet,sloupecMatice}];

(*inverzní L matice*)
Lmatice=Partition[Flatten[Lmatice1],PocCyklu+1];
inLmatice=Inverse[Lmatice];
Det[inLmatice];
inLmatice;
MatrixForm[inLmatice];


trNulovaMatice=Join[Table[{0},PocCyklu],{{1}}];
NulovaMatice=Join[Table[0,PocCyklu],{1}];

Length[NulovaMatice];
maticeIi=inLmatice.NulovaMatice;(*Matice proudů v jednoltivých lanech*)


Total[maticeIi]-maticeIi[[Length[maticeIi]]];(*Kontrola, zda je suma Ii = 1*)
Total[Delete[maticeIi,-1]];

(*Rozdělení zpět do lan*)

MaticeI=Abs[Delete[maticeIi,-1]];
MatrixForm[maticeIi];
MatrixForm[MaticeI];
Length[maticeIi];
Length[MaticeI];

Graf=Reap[For[i=1,i≤Length[MaticeI],i++,Sow[

{svazek[[i,1]],svazek[[i,2]],MaticeI[[i]]}

]];][[2,1]];

ListPointPlot3D[Graf,ColorFunction->"BlueGreenYellow",AspectRatio->Automatic,BoxRatios->Automatic,PlotStyle->PointSize[0.005]]

ListPointPlot3D[Graf,ColorFunction->"Rainbow",PlotStyle->PointSize[0.008]]

(*Výpočet Kprox*)

Suma=Reap[Sow[For[i=1,i≤PocCyklu,i++,Sow[
(((MaticeI[[i]]^2))/((Pi*(vzdalenostniMatice[[i,i]]^2))))
]
]
]][[2,1]]/.Null->Sequence[];


Kprox1=PocetLanVeSvazku*Sal*Total[Suma];
Kprox=Kprox1[[1]];
Total[Suma];

ListPlot[svazek,PlotRange->{{-1.5a-50,1.5a+50},{-0.5b-50,0.5b+50}}];
Print["Rozměry svazku: ",Style[a*100]," x ",Style[b*100]," [cm]"]
Print["kprox: ",Style[Kprox]," [%]"]

(*Tepelný model vodiče*)


ρ20=0.0384;
α20=4.03*10^(-3);
ξ20=8*10^(-7);

αs=0.5;(*Součinitel emisivity absorpce*)
ϵs=0.5;(*Součinitel emisivity*)
Qs=1120 ;(*Intenzita slunečního záření W/m2*)
σB=5.67*10^-8;(*Stefan-Boltzmanova konstanta*)

Tav=95;(*Maximalni dovolena teplota vodice*)
Ta=40;(*Teplota okoli*)
y=1000;(*Nadmořská výška*)
α=4.03*10^(-3);(*Lineární teplotní součinitel odporu*)
ξ=8*10^(-7);(*Kvadratický teplotní součinitel odporu*)


d=36.5;(*Prumer vodice [mm]*)
dp=4.12; (*Prumer dratku v posledni vrstve [mm2]*)
dfe=8.4;(*Prumer ocelove duse vodice [mm]*)
v=0.6;(*Rychlost větru*)
δ=90;(*Úhle větru působícího na vodič*)
λ=3.5; (*Dle TB601*)

ProudIn=1265; (*In [A]*)
αkonv=500;
Cp=0.241;


(*Joulovy ztraty*)
ρ=ρ20*(1+α20*(Ta-20)+ξ20*(Ta-20)^2);

Rdc=((ρ*10^-6)/Sal);
Rac=Rdc*Kprox;
Pj=Rac*ProudIn^2;

(*Sluneční záření*)
Ps=αs*Qs*d*0.001;

(*Konvekce*)

Pt=Ps+Pj;

Ts=Tav-((1/2)*(((Pt/(2*Pi*λ)))*((1/2)-(Log[(d/dfe)]*((dfe^2)/(d^2-dfe^2))))));
Tf=(1/2)*(Ts+Ta);
λf=2.368*10^-2+7.23*10^-5*Tf-2.763*10^-8*Tf^2;

γ=((1.293-(1.525*10^-4*y)+6.379*10^-9*y^2)/(1+0.00367*Tf));
μf=(17.239+(4.635*10^-2*Tf)-2.03*10^-5*Tf^2)*10^-6;
vf=((μf)/γ);
Rei=((v*d*0.001)/vf);
Rs=(dp/(2*(d-dp)));

B=If[Rs≤0.05,If[Rei≤2650,0.641,0.178],If[Rei≤2650,0.641,0.048]];(*Dle tab 4 v TB 601*)
n=If[Rs≤0.05,If[Rei≤2650,0.471,0.633],If[Rei≤2650,0.471,0.800]];(*Dle tab 4 v TB 601*)

Nu90=B*Rei^n;


Nuδ=If[δ>24,Nu90*(0.42+0.58*Sin[δ*(Pi/180)]^0.9),Nu90*(0.42+0.68*Sin[δ*(Pi/180)]^1.08)];

PcIdeal=Pi*λf*(Ts-Ta)*Nuδ;

Johrev=0.001PcIdeal*60/2;

T2=((Johrev+Cp*Ta)/Cp);

Tvys=((T2+Ta)/2);
Pc4=Pi*λf*(Ts-Tvys)*Nuδ;

Pc=If[TypSvazku==1,
If[PocetLanVeSvazku==1,PcIdeal,If[PocetLanVeSvazku==2,PcIdeal,If[PocetLanVeSvazku==3,PcIdeal,If[PocetLanVeSvazku==4,Pc4,"špatná hodnota PocetLanVeSvazku S"]]]],
If[PocetLanVeSvazku==3,PcIdeal,If[PocetLanVeSvazku==4,PcIdeal,"špatná hodnota PocetLanVeSvazku N"]]];


Pc4=Pi*λf*(Ts-Tvys)*Nuδ;

(*Radiace*)

Pr0=Pi*d*σB*ϵs*((Ts+273)^4-(Ta+273)^4)*10^-3;

(*Sálání pro 4-svazek*)
Prkoef4=(Pi-2*ArcTan[((d*0.001)/2)*a]-2*ArcTan[((d*0.001)/2)*b]-2*ArcTan[((d*0.001)/2)*(Sqrt[a^2+b^2])])*d*0.001*ϵs*σB*((Ts+273.15)^4-(Ta+273.15)^4);
(*Sáláni pro 3-svazek a plochý 4-svazek*)

Prkoef3=(Pi-4*ArcTan[((d*0.001)/2)*a])*d*0.001*ϵs*σB*((Ts+273.15)^4-(Ta+273.15)^4);

(*Sáláni pro 2-svazek*)

Prkoef2=(Pi-2*ArcTan[((d*0.001)/2)*a])*d*0.001*ϵs*σB*((Ts+273.15)^4-(Ta+273.15)^4);

(*Zatizitelnost*)


Pr=If[TypSvazku==1,
If[PocetLanVeSvazku==1,Pr0,If[PocetLanVeSvazku==2,Prkoef2,If[PocetLanVeSvazku==3,Prkoef3,If[PocetLanVeSvazku==4,Prkoef4,"špatná hodnota PocetLanVeSvazku S"]]]],
If[PocetLanVeSvazku==3,Prkoef3,If[PocetLanVeSvazku==4,Prkoef3,"špatná hodnota PocetLanVeSvazku N"]]];



Iv=(Sqrt[((Pc+Pr-Ps)/Rac)]);

p1troj=((Log[10,(a/(Ksi*((d/2)*0.001)))])/(Log[10,(1/2)]+3*Log[10,(a/(Ksi*((d/2)*0.001)))]));

p2troj=((Log[10,(1/2)]+Log[10,(a/(Ksi*((d/2)*0.001)))])/(Log[10,(1/2)]+3*Log[10,(a/(Ksi*((d/2)*0.001)))]));


p2ctyr=((Log[10,(a/(Ksi*((d/2)*0.001)))])/(4*Log[10,((2*a)/((d/2)*0.001))]+Log[10,(1/(2a))]));
p1ctyr=(1/2)(1-2*p2ctyr);

(*Trojsvazek*)
Isv3n=2*Iv+(Iv/p1troj)*p2troj;
(*print zatizitelnoist svazku je ... procent*)
pomer3=(Isv3n/(PocetLanVeSvazku*Iv));
(*Ctyrsvazek*)

Isv4n=2*Iv+2(Iv/p1ctyr)*p2ctyr;
pomer4=(Isv4n/(PocetLanVeSvazku*Iv));


IvypSv=If[TypSvazku==1,
If[PocetLanVeSvazku==1,PocetLanVeSvazku*Iv,If[PocetLanVeSvazku==2,PocetLanVeSvazku*Iv,If[PocetLanVeSvazku==3,PocetLanVeSvazku*Iv,If[PocetLanVeSvazku==4,PocetLanVeSvazku*Iv,"špatná hodnota PocetLanVeSvazku S"]]]],
If[PocetLanVeSvazku==3,Isv3n,If[PocetLanVeSvazku==4,Isv4n,"špatná hodnota PocetLanVeSvazku N"]]];




ListPlot[svazek,AspectRatio->1];
Print["Rozměry svazku: ",Style[a*100]," x ",Style[b*100]," [cm]"]

Print["kprox: ",Style[Kprox]," [%]"]

If[TypSvazku==1,
If[PocetLanVeSvazku==1,Print["Vypočtené zatížení lana: ",Style[Iv]," [A]"],
If[PocetLanVeSvazku==2,Print["Vypočtené zatížení lana: ",Style[Iv]," [A]"],
If[PocetLanVeSvazku==3,Print["Vypočtené zatížení lana: ",Style[Iv]," [A]"],
If[PocetLanVeSvazku==4,Print["Vypočtené zatížení lana: ",Style[Iv]," [A]"],
"špatná hodnota PocetLanVeSvazku S"]]]],

If[PocetLanVeSvazku==3,Print["Vypočtené zatížení krajního lana: ",Style[Iv]," [A] a vnitřního lana:" ,Style[(Iv/p1troj)*p2troj]," [A]" ],
If[PocetLanVeSvazku==4,Print["Vypočtené zatížení krajního lana: ",Style[Iv]," [A] a vnitřního lana:" ,Style[(Iv/p1ctyr)*p2ctyr]," [A]"],"špatná hodnota PocetLanVeSvazku N"]]]



Print["Jmenovité zatížení lana: ",Style[ProudIn]," [A]"]
Print["Vypočtené zatížení svazku: ",Style[IvypSv]," [A]"]