f=50;
omega=2*pi*f;
numcond=3;			%%numero dei subconduttori per fase
raggiocond=31.5/2/1000;		%%raggio dei subconduttori [m]
raggiofascio=230.94/1000;	%%raggio crf dei conduttori a fascio [m]
D_12=9;				%%distanza fasi 1-2 [metri]
D_23=9;
D_31=18;
D=(D_12*D_23*D_31)^(1/3);	%%distanza media geometrica tra le fasi
d_eq=2*(numcond*raggiocond*((raggiofascio)^(numcond-1)))^(1/numcond);

r=0;				%%resistenza chilometrica [ohm/km]
l=(0.033/numcond +0.46*log10(2*D/d_eq))*0.001; %%induttanza chilometrica [H/km]
c=24.14/(log10(2*D/d_eq))*1e-9; %% [F/km]
g=0;

z=r+i*omega*l;
y=g+i*omega*c;
%%calcolo k e Z0
k=(z*y)^0.5;
Z0=(z/y)^0.5;

Za=(2/3)*Z0;		%%linea al 150%% della pot.nat.

h=0;
for phi_a=-pi/2:pi/18:pi/2
  ++h;
  Va=380000*(cos(phi_a)+i*sin(phi_a));
  Ia=Va/Za;
  arrivo=[Va;Ia];

  for lunghezza=1:15000
    %%matrice di trasmissione
    Tpd=[cosh(k*lunghezza),Z0*sinh(k*lunghezza);
         (1/Z0)*sinh(k*lunghezza),cosh(k*lunghezza)];
    partenza=Tpd*arrivo;
    Vp(h,lunghezza)=imag(partenza(1));
    Ip(h,lunghezza)=imag(partenza(2));
  endfor

  subplot(2,1,1)
  if (h==10)		%%phi_a=0
    plot(Vp(h,1:15000),'b');
  elseif h==13		%%phi_a=30
    plot(Vp(h,1:15000),'k');
  elseif h==16		%%phi_a=60
    plot(Vp(h,1:15000),'r');
  elseif h==19		%%phi_a=90
    plot(Vp(h,1:15000),'c');
  endif
  hold on

  subplot(2,1,2)
  if (h==10)		%%phi_a=0
    plot(Ip(h,1:15000),'b');
  elseif h==13		%%phi_a=30
    plot(Ip(h,1:15000),'k');
  elseif h==16		%%phi_a=60
    plot(Ip(h,1:15000),'r');
  elseif h==19		%%phi_a=90
    plot(Ip(h,1:15000),'c');
  endif
  hold on
endfor %%ciclo for dell'argomento della tensione di arrivo

disp lol

%%Inviluppo
for lunghezza=1:15000
  inviluppoV(lunghezza)=0;
  inviluppoI(lunghezza)=0;
  for h=1:19
    if abs(Vp(h,lunghezza))>inviluppoV(lunghezza)
      inviluppoV(lunghezza)=abs(Vp(h,lunghezza));
    endif
    if abs(Ip(h,lunghezza))>inviluppoI(lunghezza)
      inviluppoI(lunghezza)=abs(Ip(h,lunghezza));
    endif
  endfor
endfor

disp plotting...

subplot(2,1,1)
plot(inviluppoV,'*');
plot(-inviluppoV,'*');
subplot(2,1,2)
plot(inviluppoI,'*');
plot(-inviluppoI,'*');
%%end Inviluppo


subplot(2,1,1)
title("Andamenti di tensione e corrente di linea aerea ideale funzionante al 150%% della propria potenza naturale")
ylabel('Tensione [V]')
axis([0 15000])
set (gca(), 'xdir', 'reverse')
hold off

subplot(2,1,2)
ylabel('Corrente [A]')
xlabel('Distanza dall''arrivo della linea [km]')
axis([0 15000])
set (gca(), 'xdir', 'reverse')
hold off

%%Salvataggio grafici
print -dpng andamenti_aerea_reale_potnat.png