{ este programa calcula los angulos de switcheo de una rama de un }
{ puente, teniendo como entrada el valor del primer armonico y los}
{ armonicos mayores que se quieren anular                         }

unit PWM_modo;

interface
uses
	xMatDefs;

var    modo:integer;
       mult:integer;
		 termino_ind:NReal;

const  MMAX = 50; { maximo valor de armonicos impares a eliminar  }
                  { en este programa es de 49                     }
       MMAX3 = 150;
		 eps:NReal=0.01;      { converge si norma del vector <= abs(eps*f1))}
       MaxInterac:integer=50; {n£mero m ximo de iteraciones antes de abortar}
		 pasominimo:NReal=1e-10;


type
	vector = array[1..MMAX] of NReal;
   fila   = array[1..MMAX] of integer;
	lvector = array[1..MMAX3] of NReal;

function fi( ni,m : integer ;var  alfa : vector) : NReal;
{    esta funcion calcula el valor de la armonica ni en funcion   }
{    del vector de angulos de switcheo alfa                       }

function calcang(var n:fila; m:integer; f1:NReal; var alfa:vector):boolean;
{retorna true si converge antes de hacer MaxInterac y el paso no se hizo }
{menor que pasominimo                                             }
{alfa: vector con los angulos                                     }
{n: vector que tiene como primer elemento 1 y luego los armonicos }
{a eliminar                                                       }
{m: numero de armonicos a eliminar + 1                            }
{f1: valor deseado de la fundamental en p.u                       }

{ calcula el contenido de armonicos dados los angulos de switcheo}
procedure conarm(var arm:lvector;var alfa:vector;m:integer;var aux:integer);
{aux: armonico hasta el que se quiere calcular el contenido}

{calcula los delta t entre conmutaciones y retorna el minimo}
function vetiemp(var dt:vector;var alfa:vector; m:integer;f:NReal):NReal;
{dt: vector de tiempos calculados}
{alfa: vectos de angulos}
{m: numero de armonicos a eliminar + 1}
{f: frecuencia de la fundamental en Hz}

procedure determinarmodo(md:integer);
{1: pwm tres estados modificado (t=0 en cero)}
{2: pwm dos estados}
{3: pwm tres estados, conmutaci¢n en t=0}

implementation

type     matriz = array[1..MMAX] of vector;

var m,i : integer;
    n : fila;
	 f1: NReal;

    alfa:vector;
    arm: lvector;                   { para poder calcular mas armonicos que}
    a:matriz;                       { los eliminados                       }
    b:vector;
    salir, otrat :boolean;

{*****************************************************************}
{    esta funcion calcula el valor de la armonica ni en funcion   }
{       del vector de angulos de switcheo alfa                    }

function fi( ni,m : integer ; var alfa : vector) : NReal;

var i,signo : integer;
	 fiaux:NReal;
begin
     signo:=-1;
     fiaux:=termino_ind;
     for i :=1 to m do
     begin
          fiaux := fiaux + mult * signo * cos( ni * alfa[i] );
          signo := -signo;
     end;
     fi := fiaux;
 end;

{***************************************************************}
{ esta procedimiento calcula el vector de funciones fi          }

procedure calcf( var f : vector; var n:fila; var alfa : vector; m : integer);
var i : integer;

begin
     for i := 1 to m do
     f[i] := fi( n[i] , m , alfa);
end;
{*****************************************************************}

Procedure sist ( var a: matriz;var x : vector; var b: vector; dim: integer);

{ resuelve el sistema lineal   a * x   +  b  =  0   }

var k,i:integer;
	 m:NReal;
{***********}
procedure inver(var A:matriz; n:integer);
(* invierte la matriz  A[n*n]  .RCH89*)
(* deter := det(A)                   *)

var
 i:matriz;
 k,p,j:integer;
 det,m:NReal;
 z:vector;

begin

for k:=1 to n do
  for j:=1 to n do
    if k=j then i[k,j]:=1
      else i[k,j]:=0;

p:=1;det:=1;
{esca1}
while p<n do
 begin
  m:=abs(a[p,p]);j:=p;
  for k:=p+1 to n do
      if abs(a[k,p])>m then
            begin
              m:=abs(a[k,p]);j:=k
            end;

  if p<>j then
          begin
           z:=a[j];a[j]:=a[p];
           a[p]:=z;z:=i[j];
           i[j]:=i[p];i[p]:=z;
           det:=-det
          end;

  if a[p,p]=0 then
               begin
                det:=0;
                p:=n
               end
      else{eliminacion}
      begin
       det:=det*a[p,p];
       for k:=p+1 to n do
         begin
          m:=-a[k,p]/a[p,p];
          for j:=p+1 to n do
            a[k,j]:=a[k,j]+m*a[p,j];
          for j:=1 to n do
            i[k,j]:=i[k,j]+m*i[p,j]
         end
       end;

  p:=p+1;


 end;
 det:=det*a[n,n];

 if det<>0 then
     begin{esca2}
      for k:=1 to n do
       begin
        for j:=1 to n do
                 i[k,j]:=i[k,j]/a[k,k];
        for j:=n downto k do
                 a[k,j]:=a[k,j]/a[k,k]
       end;
      for p:=n downto 2 do
         for k:=p-1 downto 1 do
            for j:=1 to n do
                 i[k,j]:=i[k,j]-i[p,j]*a[k,p];
       a:=i
      end;

{deter:=det}
end {inver};
{***********}
Begin
inver(a,dim);
for i:=1 to dim do
begin
     m:=0;
     for k:=1 to dim do
     m:=m+a[i,k]*(-b[k]);
     x[i]:=m;
end;
end; {sist}
{*****************************************************************}

Procedure calcgrad ( var gradiente: matriz; m: integer; var alfa: vector; var n: fila);

{este procedimiento calcula el gradiente delta_f/delta_alfa, cuyos elementos}
{ g[i,j] = - mult * (-1)^(j)* (n[i] * alfa[j]) }

var i, j, signo: integer;

Begin
for i := 1 to m do
    Begin
    signo:=-1;
    for j := 1 to m do
        Begin
        gradiente[i,j] := - mult * signo * n[i] * sin (n[i] * alfa [j]);
        signo := -signo;
        end
    end
end;
{*****************************************************************}

Function norma ( v: vector; dim: integer): NReal;

{ calcula  ((v[1])^2 + ... + (v[dim])^2)^ (1\2) }

var i: integer;
	 normaux :NReal;

Begin
normaux := 0;
for i := 1 to dim do
    Begin
    normaux := normaux + v[i] * v[i]
    end;
norma := sqrt( normaux )
end;  { norma }

{*****************************************************************}

function calcang( var n : fila; m: integer; f1: NReal; var alfa:vector):boolean;

{ esta funcion calcula los angulos de switcheo con los siguientes datos:}
{n: vector que tiene como primer elemento 1 y luego los armonicos }
{   a eliminar                                                    }
{m: numero de armonicos a eliminar + 1                            }
{f1: valor deseado de la fundamental en p.u                       }
{modo: 1 2 o 3.}

var f0, alfa0, dalfa, alfaux: vector;
    grad :matriz;
    contit, i : integer;
    ord, fin, converge{, tecla} : boolean;
	 ant,nact,nant,paso : NReal;

Begin
for i := 1 to m do alfa0[i] := Pi * i/2 / (m + 1);

contit := 0;
fin := FALSE;
converge:=true;
calcf( f0, n, alfa0, m);
f0[1] := f0[1] - f1;
nant:=norma(f0,m);

while ( not fin and converge) do
    Begin
      if ( nant <= abs(eps*f1)) then fin := TRUE
      else
        Begin
           paso:=1;
           calcgrad ( grad, m, alfa0, n);
           sist ( grad, dalfa, f0, m);
           repeat
                for i := 1 to m do alfaux[i] := alfa0[i] + paso * dalfa[i];
                calcf(f0,n,alfaux,m);
					 f0[1] := f0[1]  - f1;
                paso := paso/2;
                i := 1;
                ant := 0;
                ord := TRUE;
                while ( ord and (i < m + 1)) do
                Begin
                      ord := ( alfaux[i] > ant);
                      ant := alfaux[i];
                      i := i +1;
                end;
                ord := ord and (alfaux[m] < pi/2);
                nact := norma(f0,m);

			  until(((nact < nant) and ord) or (paso < pasominimo));

           for i:=1 to m do alfa0[i] := alfaux[i];
           nant := nact;
        end; {else}
      contit := contit + 1;
      if contit=MaxInterac then converge:=false;
    end; {while}
    for i := 1 to m do alfa[i] := alfa0[i];
    calcang:=converge;
end; { calcang }

{*****************************************************************}

{ calcula el contenido de armonicos dados los angulos de switcheo}
procedure conarm(var arm:lvector;var alfa:vector;m:integer;var aux:integer);

var i:integer;

begin
     aux:=2*(aux div 2) +1; {para que aux quede impar}
     for i:=0 to aux do
     begin
    //rch@080328, OJO esto no estaba comentado
    //  i:=i+1;
          arm[i] := fi(i, m, alfa)/i;
     end;
end;
{*************************************************************************************}

{calcula los delta t entre conmutaciones y retorna el minimo}
function vetiemp(var dt:vector;var alfa:vector; m:integer;f:NReal):NReal;
var i,j,k:integer;
    auxc:char;
    t :vector;
    min:NReal;

Begin
     for i:=1 to m do t[i] := alfa[i] /f/(2*pi);
     t[m+1] := (1/(2*f)) - t[m];
     if modo=1 then dt[1]:=2*t[1] else dt[1]:=t[1];
     min:=dt[1];
     for i:=2 to m do
     begin
          dt[i]:= t[i] - t[i-1];
          if dt[i]<min then min:=dt[i];
     end;
     dt[m+1] := t[m+1] - t[m];
     if dt[m+1]<min then min:=dt[m+1];
     vetiemp:=min;
end;

procedure determinarmodo(md:integer);
{1: pwm tres estados modificado (t=0 en cero)}
{2: pwm dos estados}
{3: pwm tres estados, conmutaci¢n en t=0}
begin
     modo:=md;
     case modo of
     1: begin termino_ind:=0;mult:=-1; end;
     2: begin termino_ind:=1;mult:=2; end;
     3: begin termino_ind:=1;mult:=1; end;
     end;
end;

begin
determinarmodo(2); {por si se olvidan es el default}
end.
