unit uVhorarioMet;
interface
uses
  sysutils, xmatdefs, MatReal,
  Classes, Graphics, math01;

// lector de velocidad de viento horaria en formato de Metereología.


const
  HUECO_DE_DATO= -11111;  // HUECO_DE_DATO indica que no hay medida
  FUERA_DE_RANGO= -11112; // indica que estamos fuera del horizonte de tiempo medido
  UMBRAL_DATO_VALIDO= 500;

type

TProc_CalcXY = procedure ( var X, Y : integer; Latitud, Longitud: NReal ) of object;



TArchiVhoraria = class
  archi: string;

  // componentes de la velocidad en m/s
	velsx, velsy: TDAOfNReal;

  dtInicial: TDateTime;
  dtFinal: TDateTime; // de inicio de la hora siguiente a la última
  cnt_999: integer;
  latitud, longitud: NReal;
  iActual: integer;
  iFinal: integer; // i de la hora siguiente a la última
  inicializada: boolean;
  xbof: boolean; // before of file
  xeof: boolean; // end of file

  constructor Create( archi: string; latitud, longitud: NReal);

  procedure Free; virtual;
  procedure init_iActual( dtAhora: TDateTime );
  procedure inc_iActual;
  procedure Draw(
    c: TCanvas; calcxy: TProc_CalcXY;
    modo: integer;
    pixByms: double // factor de conversión [Pixeles/(m/s)]
     );

  function vx: NReal; // retorna la componente vx para la iActual
  function vy: NReal; // retorna la componente vy para la iActual;

// Retorna el factor de planta para un Vestas80-2MW
// Se supone un único aero-generador por lo que no hay interferencia
// entre los aerogeneradores.
// Retorna la potencia media de las horas en las que hay datos dividida
// por 2MW.
  function FactorPlanta_V80_2MW( speedUP: NReal; vMaxParaSalida: NReal ): NReal;

// Calcula el SpeedUp necesario para llegar al factor de planta dado (fp)
// considerando la curva de potencia del Vestas V80.
  function SpeedUPparaFP( fp: NReal ): NReal;

// Multiplica la serie por el valor fp.
// no modifica los valores HUECOS
  procedure multiplicarPorReal( fp: NReal );
private
  procedure chequear_xbeof;
end;


TVientosHorarios= class
  vientos: TList;
  camino: string;
  cnt_xbof: integer; // before of file
  cnt_xeof: integer; // end of file

  constructor Create;
  procedure Free; virtual;
  procedure init_iActual( dtAhora: TDateTime );
  procedure inc_iActual;
  procedure Draw(
    c: TCanvas; calcxy: TProc_CalcXY;
    modo: integer;
    pixByms: double // factor de conversión [Pixeles/(m/s)]
     );

  // retorna el mayor de los dtInicial de cada archivo.
  function dtInicial: TDateTime;

  // retorna el menor de los dtFinal de cadaarchivo.
  function dtFinal: TDateTime;

  // si el archivo de datos no comienza "\" o con "?:\" le pegamos
  // al inicio camino+'\' a los archivos pasados en AddArchiDatos
  procedure setCaminoADatos( camino: string );
  procedure AddArchiDatos( archi: string; latitud, longitud: NReal );

  procedure WriteSeriesToXLT( archi: string );
end;



// potencia en kW del Vestas V80-2MW
function PotV80_kW( v: NReal; speedUP: NReal; vMaxParaSalida: NReal ): NReal;

// borra desde el inicio hasta la última / y luego desde el final hasta el punto.
function getNombreArchi( s: string ): string;

implementation

{$IFDEF dbg_dir2grad}
var
  fsal: textfile;
{$ENDIF}

var
  vientoActivoParaDicot: TArchiVHoraria;
  fpParaDicot: NReal;
  curvaV80_Pot, curvaV80_Vel: TVectR;


function funcParaDicot( x: NReal ): NReal;
var
  fp: NReal;
begin

  fp:=
    vientoActivoParaDicot.FactorPlanta_V80_2MW(
        x,
        1000000); // este valor es para que no saque el molino por sobrevelocidad.

  result:= fp - fpParaDicot;
end;





// potencia del Vestas V80-2MW
function PotV80_kW( v: NReal; speedUP: NReal; vMaxParaSalida: NReal   ): NReal;
var
  kr: NReal;
  vel, P: NReal;
begin
  vel:= v* speedUP;

  if ( vel > vMaxParaSalida ) then
  begin
    result:= 0;
    exit;
  end;

  kr:= curvaV80_Vel.inv_interpol( vel );
  P:= curvaV80_Pot.interpol( kr );
  result:= P;
end;

function PotV80_mps( P: NReal  ): NReal;
var
  kr: NReal;
  vel: NReal;
begin
  if ( P <= 0 ) then
  begin
    result:= 0;
    exit;
  end;

  kr:= curvaV80_Pot.inv_interpol( P );
  vel:= curvaV80_Vel.interpol( kr );
  result:= vel;
end;



// borra desde el inicio hasta la última \ y luego desde el final hasta el punto.
function getNombreArchi( s: string ): string;
var
  i: integer;
  ts: string;
begin
  ts:= s;
  while pos('\', ts ) > 0 do
    delete( ts, 1, pos( '\', ts ) );
  i:= pos( '.', ts );
  if i > 0 then
    delete( ts, i, length( ts ) - i + 1 );
  result:= ts;
end;

function Dir1ToGrad( c: char ): double;
begin
  case c of
    'N': result:= 270;
    'E': result:= 180;
    'S': result:= 90;
    'W': result:= 0;
  else
    result:= 0;
  end;
end;

function Dir2ToGrad( c1, c2: char ): double;
var
  ang1, ang2: double;
begin
  ang1:= Dir1ToGrad( c1 );
  ang2:= Dir1ToGrad( c2 );
  if abs( ang1 - ang2 ) > 91 then
    result:= (ang1+ang2+360)/2.0
  else
    result:= (ang1+ang2)/2.0;
end;


function Dir3ToGrad( c1, c2, c3: char ): double;
var
  ang1, ang2: double;
begin
  ang1:= Dir1ToGrad( c1 );
  ang2:= Dir2ToGrad( c2, c3 );
  if abs( ang1 - ang2 ) > 91 then
    result:= (ang1+ang2+360)/2.0
  else
    result:= (ang1+ang2)/2.0;
end;

function DirToGrad( dir: string ): double;
var
  n: integer;
  res: double;
begin
  n:= length( dir );
  case n of
    0: res:= 0;
    1: res:= Dir1ToGrad( dir[1] );
    2: res:= DIr2ToGrad( dir[1], dir[2] );
    3: res:= Dir3ToGrad( dir[1], dir[2], dir[3] );
  else
    raise Exception.Create('DirToGrad, largo dir inválido n:'+IntToStr( n ));
  end;
  result:= res;

{$IFDEF dbg_dir2grad}
  writeln( fsal, dir, #9, res );
{$ENDIF}

end;

function DirToRad( dir: string ): double;
var
  ang: double;
begin
  ang:= DirToGrad( dir );
  if ang > 180 then
    ang:= ang -360;
  result:= ang / 180.0* pi;
end;

procedure TArchiVHoraria.chequear_xbeof;
begin
  if iActual >= iFinal then
    xeof:= true
  else
    xeof:= false;

  if iActual < 0 then
    xbof:= true
  else
    xbof:= false;
end;

procedure TArchiVhoraria.init_iActual( dtAhora: TDateTime );
begin
  iActual:= trunc( (( dtAhora + 10.0/(24.0*60.0) ) - dtInicial)*24.0 );
  chequear_xbeof;
  inicializada:= true;
end;

procedure TArchiVhoraria.inc_iActual;
begin
  inc( iActual );
  chequear_xbeof;
end;


constructor TArchiVhoraria.Create( archi: string; latitud, longitud: NReal );
var
	f: textfile;
	r: string;
	Anio, Mes, Dia, Hora: integer;
	Dir: string;
	vel: double;
	cnt_lineas: integer;
	buscando: boolean;
	k: integer;
	DiaAnt: integer;
  timeStamp: string;
  ang: double;
  dirant: string;
begin
  dirant:= 'S';

writeln('Procesando archivo: ', archi );

  inherited Create;
  inicializada:= false;

  self.archi:= archi;
  self.latitud:= latitud;
  self.longitud:= longitud;
	vel:= 0;

  assignFile( f, archi );
  filemode:= 0;

  {$I-}
  reset(f );
  {$I+}
  if ioresult <> 0 then
    raise Exception.Create( 'No encontré el archivo: '+archi );

(* PRIMERO CONTAMOS LAS LINEAS DEL ARCHIVO *)
  cnt_lineas:= 0;
  cnt_999:= 0;
  buscando:= true;
  while not eof( f ) and buscando do
  begin
    readln( f, r );
    if (length( r ) = 14) then
    begin
      if pos('999', r )>0 then
      begin
        inc( cnt_999 );
      end;
      inc( cnt_lineas )
    end
    else
    begin
      writeln('Terminé por encontrar: ['+r+']' );
      buscando:= true;
    end;
  end;
  reset(f );

  DiaAnt:= -1;
  setlength( velsx, cnt_lineas );
  setlength( velsy, cnt_lineas );

  for k:= 0 to cnt_lineas-1 do
  begin
    readln( f, r );

//writeln( k: 5,'>', r );
    if length( r ) > 8 then
    begin
      timeStamp:= copy( r, 1, 8 );
      Anio:= 1900+ StrToInt( copy( r, 1, 2 ));
      Mes:= StrToInt( copy( r, 3, 2 ));
      Dia:= StrToInt( copy( r, 5, 2 ));
      Hora:= StrToInt( copy( r, 7, 2 ));
      Dir:= copy( r, 9, 3 );

      if k= 0 then
        dtInicial:= encodedate( Anio, Mes, Dia )+ Hora/24.0;

      vel:= StrToInt( copy( r, 12, 3));

      if vel > 990 then
      begin
        velsx[k]:= -HUECO_DE_DATO; // lo marco así para filtrarlo al analizar
        velsy[k]:= -HUECO_DE_DATO; // lo marco así para filtrarlo al analizar
      end
      else
      begin
        dir:= trim( dir );
        if (dir= '') and ( vel <> 0 ) then
        begin
          writeln('Ojo, vel: '+FloatToStr( vel )+' y dir: vacío. ' );
          writeln( 'archi: ', archi );
          writeln('klin: ', k );
          // asumo ultima dirección.
          dir:= dirant;
        end;
        dirant:= dir;

        if dir = '' then
          ang:= 0
        else
          ang:= DirToRad( trim(Dir) );
        vel:= vel / 3.6; // lo pasamos de km/h a m/s
        velsx[k]:= vel * cos( ang );
        velsy[k]:= vel * sin( ang );
      end;
    end;
  end;
  closefile( f );

   dtFinal:= dtInicial + ( 1.0/ 24.0* cnt_lineas );
   iFinal:= cnt_lineas;
end;




procedure TArchiVhoraria.multiplicarPorReal( fp: NReal );
var
  k: integer;
begin
  for k:= 0 to high( velsx ) do
    if ( abs(velsx[k]) < UMBRAL_DATO_VALIDO ) and ( abs(velsy[k]) < UMBRAL_DATO_VALIDO) then
    begin
      velsx[k]:= velsx[k]* fp;
      velsy[k]:= velsy[k]* fp;
    end;
end;


function TArchiVhoraria.FactorPlanta_V80_2MW( speedUP: NReal; vMaxParaSalida: NReal ): NReal;
var
  energia: NReal;
  cnt: integer;
  v, P: NReal;
  k: integer;
begin
  energia:= 0;
  cnt:= 0;
  for k:= 0 to high( velsx ) do
    if ( abs(velsx[k]) < UMBRAL_DATO_VALIDO ) and ( abs(velsy[k]) < UMBRAL_DATO_VALIDO) then
    begin
      v:= sqrt( sqr( velsx[k] )+ sqr( velsy[k] ) );
      P:= PotV80_kW( v, speedUP, vMaxParaSalida )/ 1000;
      Energia:= Energia + P;
      inc( cnt );
    end;

  if cnt > 0 then
    result:= energia/(cnt * 2.0 )
  else
    result:= 0;

end;

// Calcula el SpeedUp necesario para llegar al factor de planta dado (fp)
// considerando la curva de potencia del Vestas V80.
function TArchiVhoraria.SpeedUPparaFP( fp: NReal ): NReal;
var
  Root, fAtRoot: NReal;
  NoOfIts: word;
  converged: boolean;
begin
  fpParaDicot:= fp;

  vientoActivoParaDicot:= Self;
  Dicot(
  	funcParaDicot,							{funci¢n a anular}
	  0, 10, 0.001,          {extremos y tolerancia}
	  1000,					{n£mero m ximo de iteraciones}
	  Root,fAtRoot,     {ra¡z y f(ra¡z)}
	  NoOfIts,          {n£mero de iteraciones realizadas}
	  converged);		{valid‚s del resultado}

  if not converged then
    raise Exception.Create('No logré calcular SpeedUP para fp: '+FloatToStrF( fp, ffFixed, 4, 2 ) );
  result:= root;
end;



procedure TArchiVhoraria.Free;
begin
  setlength( velsx, 0 );
  setlength( velsy, 0 );
  inherited Free;
end;

procedure TArchiVHoraria.Draw(
  c: TCanvas;
  calcxy: TProc_CalcXY;
  modo: integer;
  pixByms: double );
var
  x, y: integer;
  vx, vy: integer;
begin
  calcxy(x, y, latitud, longitud );

  if modo= 0 then
  begin
    c.Brush.Color:= clRed;
    c.FillRect( Rect(x-10, y-10, x+10, y+10 ));
  end;

  if modo= 2 then
    c.Pen.Color:= clOlive
  else
    c.Pen.Color:= clYellow;
  c.Pen.Width:= 3;

  if ( xeof or xbof ) then exit;


  if ( abs(velsx[iActual ]) ) < 500 then
  begin
    vx:= trunc( pixByms * velsx[ iActual ] + 0.5 );
    vy:= trunc( pixByms * velsy[ iActual ] + 0.5 );
    c.MoveTo( x, y );
    c.LineTo( x+vx, y-vy);
  end;
 end;


function TArchiVHoraria.vx: NReal; // retorna la componente vx para la iActual
begin
  if xbof or xeof then
    result:= FUERA_DE_RANGO
  else
    result:= velsx[ iActual ];
end;

function TArchiVHoraria.vy: NReal; // retorna la componente vy para la iActual;
begin
  if xbof or xeof then
    result:= FUERA_DE_RANGO
  else
    result:= velsy[ iActual ];
end;



(********************************************************
Métodos de TVientosHorarios
=========================================================*)

constructor TVientosHorarios.Create;
begin
  inherited Create;
  vientos:= TList.Create;
  camino:= '';
end;

procedure TVientosHorarios.Free;
var
  k: integer;
  av: TArchiVHoraria;
begin
  for k:= 0 to Vientos.Count - 1 do
  begin
    av:= vientos[k];
    av.Free;
  end;
  vientos.Free;
  inherited Free;
end;

procedure TVientosHorarios.init_iActual( dtAhora: TDateTime );
var
  k: integer;
  av: TArchiVHoraria;
begin
  cnt_xbof:= 0;
  cnt_xeof:= 0;
  for k:= 0 to Vientos.Count - 1 do
  begin
    av:= vientos[k];
    av.init_iActual( dtAhora );
    if av.xbof then inc( cnt_xbof );
    if av.xeof then inc( cnt_xeof );
  end;
end;

procedure TVientosHorarios.inc_iActual;
var
  k: integer;
  av: TArchiVHoraria;
begin
  cnt_xbof:= 0;
  cnt_xeof:= 0;
  for k:= 0 to Vientos.Count - 1 do
  begin
    av:= vientos[k];
    av.inc_iActual;
    if av.xbof then inc( cnt_xbof );
    if av.xeof then inc( cnt_xeof );
  end;
end;
                                    // mapa grande>-
procedure TVientosHorarios.Draw(
  c: TCanvas;
  calcxy: TProc_CalcXY;
  modo: integer;
  pixByms: double );
  
var
  k: integer;
  av: TArchiVHoraria;
begin
  for k:= 0 to Vientos.Count - 1 do
  begin
    av:= vientos[k];
    av.Draw( c, calcxy, modo, pixByms );
  end;
end;

// si el archivo de datos no comienza "\" o con "?:\" le pegamos
// al inicio camino+'\' a los archivos pasados en AddArchiDatos
procedure TVientosHorarios.setCaminoADatos( camino: string );
begin
  self.camino:= camino;
end;

procedure TVientosHorarios.AddArchiDatos( archi: string; latitud, longitud: NReal );
var
  k: integer;
  archix: string;
  av: TArchiVHoraria;
begin
  k:= pos( '\', archi );
  if ( k= 1) or ( (k= 3) and (archi[2]=':')) then
    archix:= archi
  else
    archix:= camino+'\'+archi;

  av:= TArchiVHoraria.Create( archix, latitud, longitud );
  vientos.Add(av);
end;

// retorna el mayor de los dtInicial de cada archivo.
function TVientosHorarios.dtInicial: TDateTime;
var
  k: integer;
  av: TArchiVHoraria;
  res: TDateTime;
begin
  if Vientos.Count = 0 then
  begin
    result:= 0;
    exit;
  end;

  av:= vientos[0];
  res:= av.dtInicial;
  for k:= 1 to Vientos.Count - 1 do
  begin
    av:= vientos[k];
    if  av.dtInicial > res then
      res:= av.dtInicial;
  end;
  result:= res;
end;

// retorna el menor de los dtFinal de cadaarchivo.
function TVientosHorarios.dtFinal: TDateTime;
var
  k: integer;
  av: TArchiVHoraria;
  res: TDateTime;
begin
  if Vientos.Count = 0 then
  begin
    result:= 0;
    exit;
  end;

  av:= vientos[0];
  res:= av.dtFinal;
  for k:= 1 to Vientos.Count - 1 do
  begin
    av:= vientos[k];
    if  av.dtFinal < res then
      res:= av.dtFinal;
  end;
  result:= res;
end;


procedure TVientosHorarios.WriteSeriesToXLT( archi: string );
var
  fsal: textfile;
  dt1, dt2, dtActual: TDateTime;
  kHora: integer;
  iViento: integer;
  viento: TArchiVhoraria;
  fp{, varp}: NReal;
  speedup: array of NReal;
  vx, vy: NReal;
  P, v: NReal;
  cnt: integer;
  acumP: NReal;

begin
  assignfile( fsal, archi );
  rewrite( fsal );

  dt1:= dtInicial;
  dt2:= dtFinal;

  if dt1 >= dt2 then
  begin
    writeln(fsal, 'dt1: ', DateTimeToStr( dt1 ) );
    writeln(fsal, 'dt2: ', DateTimeToStr( dt2 ) );
    writeln(fsal, 'NO hay una venta común a todos los archivos. No es posible simular' );
    closefile( fsal );
    exit;
  end;

// los ponemos todos al inicio.
  init_iActual( dt1 );
  kHora:= 0;
  dtActual:= dt1 + kHora / 24.0;


// escribimos el encabezado
  write( fsal, '---' );
  for iViento:= 0 to vientos.count-1 do
  begin
    viento:= vientos[iViento];
    write( fsal, #9, getNombreArchi( viento.archi ) );
  end;
  writeln( fsal );

// escribimos los factores de planta y speedUP
  write( fsal, 'fp' );
  setlength( speedup, vientos.Count );
  for iViento:= 0 to vientos.count-1 do
  begin
    viento:= vientos[iViento];
    fp:= viento.FactorPlanta_V80_2MW( 1.0, 25 );
    write( fsal, #9, fp:4:2 );
  end;
  writeln( fsal );


// escribimos los speedUP para FP=0.3
  write( fsal, 'SU03' );
  setlength( speedup, vientos.Count );
  for iViento:= 0 to vientos.count-1 do
  begin
    viento:= vientos[iViento];
    speedUP[iViento]:=  viento.SpeedUPparaFP( 0.3 );
    write( fsal, #9, speedUP[iViento]:4:2 );
  end;
  writeln( fsal );

// escribimos la Latitud
  write( fsal, 'Latitud [º]' );
  for iViento:= 0 to vientos.count-1 do
  begin
    viento:= vientos[iViento];
    write( fsal, #9, viento.latitud:4:1 );
  end;
  writeln( fsal );

// escribimos la Longitud
  write( fsal, 'Longitud [º]' );
  for iViento:= 0 to vientos.count-1 do
  begin
    viento:= vientos[iViento];
    write( fsal, #9, viento.longitud:4:1 );
  end;
  writeln( fsal );

  writeln( fsal, 'Potencias  fp03' );
// escribimos el encabezado 3
  while cnt_xeof = 0 do
  begin
    dtActual:= dt1 + kHora / 24.0;
    write( fsal, DateTimeToStr( dtActual  ) );


    cnt:= 0;
    acumP:= 0;
    for iViento:= 0 to vientos.count-1 do
    begin
      viento:= vientos[iViento];
      vx:= viento.vx;
      vy:= viento.vy;
      if ( abs( vx ) > 1000 ) or ( abs( vy ) > 1000 ) then
        P:= HUECO_DE_DATO
      else
      begin
        v:= sqrt( sqr( vx ) + sqr( vy ) );
        P:= PotV80_kW( v, speedUP[iViento], 2500);
        acumP:= acumP + P;
        inc( cnt );
      end;
      write( fsal, #9, P:8:3);
    end;
    if cnt > 0 then
      acumP:= acumP / cnt;

    v:= PotV80_mps( acumP );
    writeln( fsal, #9, v: 8:2 );
    inc_iActual;
    inc( kHora );
  end;

  closefile( fsal );
end;

initialization
{$IFDEF dbg_dir2grad}

assign( fsal, 'x.xlt' );
rewrite( fsal );
{$ENDIF}

curvaV80_Pot:= TVectR.Create_Init( 7 );
curvaV80_Vel:= TVectR.Create_Init( 7 );

with curvaV80_Pot do
begin
  pv[1]:= 0;
  pv[2]:= 3.3;
  pv[3]:= 200;
  pv[4]:= 800;
  pv[5]:= 1400;
  pv[6]:= 2000;
  pv[7]:= 2001;
end;

with curvaV80_vel do
begin
  pv[1]:= 0;
  pv[2]:= 3.3;
  pv[3]:= 5.4;
  pv[4]:= 8.4;
  pv[5]:= 10.4;
  pv[6]:= 13.8;
  pv[7]:= 113.0;
end;

finalization
{$IFDEF dbg_dir2grad}
closefile( fsal );
{$ENDIF}

curvaV80_pot.Free;
curvaV80_vel.Free;
end.
