{+doc
+NOMBRE: it4
+CREACION: febrero 93
+AUTORES: BZ
+REGISTRO:
+TIPO: Unidad Pascal.
+PROPOSITO: Transformada de BZ
+PROYECTO: SOLAR.
+DESCRIPCION: Cambio de it2 para que de datos horarios.
	Antes de utilizar los datos de un d¡a hay que llamar la funcion
	TransformadaAPlanoInclinado para que se carge el vector EnergiaHora con
	los valores correspondientes.

-doc}

{ utiliza el archivo de datos ADATOS }

UNIT IT4;

{$MODE Delphi}

INTERFACE
uses
  SysUtils, dialogs;

type
	dia=1..366;


type
	PVectDatHr=^TVectDatHr;
	TVectDatHr=  array[0..23] of real;

var
	EnergiaHoraPtr: PVectDatHr;

function NuevoDia(
	ndia:dia;   { D¡a del a¤o. 1 = primero de enero }
	beta:real;       { Inclinaci¢n del panel [radianes] }
	modelo:string; {modelo de rad. solar, existen 3 modelos
						 'Liu-Jordan' 'Perez' y 'HDKR'}
	Hg: real    { Energ¡a global en el d¡a (ndia) sobre
               plano horizontal en la superficie de la
               tierra. [WH/m2] }

	 ):real; { [WH/m2], energía diaria en el plano inclinado }



function NuevoDia_HgCalcm2(
	ndia:dia;   { Dia del anio. 1 = primero de enero }
	beta:real;       { Inclinacion del panel [radianes] }
	modelo:string; {modelo de rad. solar, existen 3 modelos
						 'Liu-Jordan' 'Perez' y 'HDKR'}
	Hg: real { [cal/cm2] }
	 ):real; { [Wh/m2] }

function Cal_A_KWH( cal:real):real;
function KWH_A_Cal( KWH:real):real;

procedure IT_Init( beta:real);


procedure logwrite( s: string );

IMPLEMENTATION


const
  pi =     3.141592;   {n£mero pi}
  latitud  =    -34.9*pi/180; {latitud de Montevideo}
  albedo=  0.25; { albedo para el pasto }
  step=    0.02;
  Isc=     1367;      { Valor medio de la intensidad solar
                      extraterrestre plano normal [W/m2] }
  arckt= 'pasb.kt';



type

	TMatIndices= array[1..390, 1..30]  of integer;
	PMatIndices= ^TMatIndices;

  meses= 1..12;
  vector= packed array[1..14] of real;

  vector2= array[dia] of real;


  registrodedatos  = packed record
                    ktdato : real ;
                    fraccionglobal : vector ;
                    fracciondirecta :vector ;
                    end;
  datos = file of registrodedatos ;



var
   registro : registrodedatos;
   archivo  : datos ; // se abre al  inicio.
   dias:dia;
   mes:meses ;

   Isn,delta,Ws,Ho,
   sin_Ws,cos_delta,sin_delta :vector2 ;
   total: array[1..390] of integer;
   ref: array[1..12] of integer;
	 indicePtr: PMatIndices;
   efe: array[1..8,1..6] of real;

	 tii,tff,
   grupokt:integer;

   sin_L,cos_L,cos_lbeta,sin_lbeta,cos_beta,sin_beta,sin_betac,
	 ti,tf,wii,wff,kt, I0hora: real;

	 model:integer; {modelo para calcular la
							 radiacion en plano inclinado
							 1-> Liu-Jordan
							 2-> HDKR
							 3-> Perez}
{--------------------------------------------}

Function CALporCM2_A_KJporM2( x:real):real;
begin
   CALporCM2_A_KJporM2:=x*41.86;
end;
{--------------------------------------------}

Function Tan(x:real):real;
begin
  Tan:= sin(x)/cos(x)
end;
{--------------------------------------------}

Function Acos(x:real):real;
begin
  Acos:= -Arctan(x/sqrt(1-x*x)) +pi/2 ;
end;
{--------------------------------------------}

Function Cal_A_KWH( cal:real):real;
begin
   Cal_A_KWH:=cal*1.162E-6
end;

Function KWH_A_Cal( KWH:real):real;
begin
	 KWH_A_Cal:=KWH/1.162E-6
end;

{--------------------------------------------}

Function Ktf(Hg:real):real;
begin
	 Ktf:=Hg/Ho[dias];
end;
{--------------------------------------------}

Procedure Busqueda(Hg:real;var registro:registrodedatos);
var
    grupoktmedio,elegido,fila: integer;
begin

{} {append(arckts);}

	 {grupoktmedio surge de ktmax=1 de donde el indicePtr^ del grupo con kt=1/2
    sera (1/2/step +0.5)}


 grupoktmedio:= trunc ((1/2/step +0.5)/2 ) +30*(mes-1) ;
 kt:= Ktf(Hg);
 grupokt:= trunc(kt/2/step+1) +30*(mes-1) ;

 {el grupokt puede no contener filas asignadas}
    while total[grupokt] = 0 do
     begin
{}  {   writeln(arckts,dias,'   ',grupokt,'   ',total[grupokt],'   ',kt:1:3 );}
       if grupokt > grupoktmedio then
     grupokt:=grupokt-1
       else
         grupokt:= grupokt+1
       end;

    {busca un numero entre 1 y el total del grupokt}
    elegido:= trunc(random*total[grupokt]) + 1 ;
		fila:= indicePtr^[ grupokt,elegido ];
    seek(archivo,fila-1);
    read(archivo,registro);

{} { writeln(arckts,dias,'   ',grupokt,'   ',fila,'   ',
registro.ktdato:1:3,'   ',kt:1:3 );}
{}{  writeln(dias,'   ',grupokt,'   ',fila,'   ',
registro.ktdato:1:3,'   ',kt:1:3 );}

{}{ close(arckts);}
end;
{--------------------------------------------}

Procedure Salida_Puesta(var ti,tf,wii,wff:real;var tii,tff:integer);
var
   td:real;
begin
  td:= 2*180/pi/15 *ws[dias];
  ti:= 12- td/2;
  tii:= trunc(ti);
  tf:= 12+td/2;
  tff:= trunc(tf+1);
  wii:= pi/12*(trunc(ti)+ 1 -12);
  wff:= pi/12*(trunc(tf)-12) ;
end;
{--------------------------------------------}

Procedure Whora(hora:integer;var wmed,wdif:real);

var
w1,w2 :real ;

begin
   if hora=tii then
        begin
		  w1:=-ws[dias];
        w2:= wii;
        end
   else
      if hora= trunc(tf) then
         begin
         w1:=wff;
			w2:=ws[dias];
         end
      else
         begin
         w1:=wii+pi/12*(hora-tii-1);
         w2:= w1+pi/12
         end;
  wmed:= (w1+w2)/2;
  wdif:= w2-w1;

end;
{--------------------------------------------}

Procedure Verificacion(iDia: Dia; hora,wmed,wdif:real;var global,directa,difusa:real);
var
    kthora: real;

begin
  I0hora:= 3.6/41.86 * Isn[iDia]*( cos_l*cos_delta[iDia]*2*sin(wdif/2)*cos(wmed)
            + wdif* sin_l*sin_delta[iDia] ) ; { cal/cm2}
  kthora:= global/I0hora;
  if kthora>1 then
   begin


   directa:= directa/global* 0.9*I0hora;
    global:= 0.9*I0hora ;
   difusa:= global- directa;
   end;
end;
{--------------------------------------------}

Function Liu_Jordan(cosinc,cosho,global,directa,difusa:real):real;

var
  reflejada:real;
begin

{ if directa<0 then writeln('directa <0 ');
 if cosho<0 then writeln('cosho <0  ');}

 if cosinc <0 then
  directa:=0
 else
  directa:= cosinc/cosho*directa;

  difusa:= (1+cos_beta)/2*difusa;
  reflejada:= albedo*(1-cos_beta)/2*global;

  Liu_Jordan:= directa + difusa +reflejada;
{  writeln(trunc(directa),'  ',trunc(difusa),'    ',trunc(reflejada),'   ');}
end;
{--------------------------------------------}

Function Hdkr(cosinc,cosho,global,directa,difusa:real):real;
{HAY DAVIES KLUCHER REINDL}

var
f,ai,reflejada : real;

begin
 if global=0 then f:=1 else f:= sqrt(directa/global);
 ai:= directa/I0hora;

 if cosinc <0 then
  directa:=0
 else
  directa:= cosinc/cosho*directa;

  difusa:= (ai*cosinc/cosho + (1-ai)*(1+cos_beta)/2 *(1+f*sin_betac))*difusa;
  reflejada:= albedo*(1-cos_beta)/2*global;
  HDKR:= directa + difusa +reflejada;
end;

{--------------------------------------------}

Function Perez(beta,wdif,cosinc,cosho,
                global,directa,difusa:real):real;

var

aa,bb,ab,exn,del,directan,titaz,eps,f1,f2,tit3  : real;
reflejada: real;
ij: integer;

begin
if cosinc > 0 then aa:= cosinc else aa:= 0;
if cosho > cos(pi*85/180) then bb:= cosho
else bb:= cos(pi*85/180);

ab:=aa/bb ;
if beta=0 then ab:=1 ;

exn:= Isn[dias]*wdif;
del:= difusa/(exn * 3.6 * cosinc);
directan:= directa/cosho;
titaz:= Acos(cosho) ;
tit3:=titaz*titaz*titaz;

if difusa > 1.0E-8 then
begin
  eps:= ((difusa+directan)/difusa + 1.04108*tit3 )/(1+1.04108*tit3);
  {para difusa nula hor.dara una componente nula de difusa en el plano inclinado}

  if (eps <= 1.065) then ij:= 1
  else if (eps > 1.065) and (eps <= 1.230) then ij:= 2
  else if (eps > 1.230) and (eps <= 1.500) then ij:= 3
  else if (eps > 1.500) and (eps <= 1.950) then ij:= 4
  else if (eps > 1.950) and (eps <= 2.800) then ij:= 5
  else if (eps > 2.800) and (eps <= 4.500) then ij:= 6
  else if (eps > 4.500) and (eps <= 6.200) then ij:= 7
  else ij:= 8;

  {} {writeln(eps);}

  f1:= efe[ij,1] + efe[ij,2]*del + titaz*efe[ij,3];
  if f1<0 then f1:=0;
  f2:= efe[ij,4] + efe[ij,5]*del + titaz*efe[ij,6];
  difusa:= ( (1-f1)*(1+cos_beta)/2 + f1*ab + f2*sin_beta )*difusa;
end
else
  difusa:= 0;

 if cosinc <0 then
   directa:=0
 else
   directa:= cosinc/cosho*directa;

reflejada:= albedo*(1-cos_beta)/2*global;

Perez:=  directa +difusa +reflejada ;
end;
{--------------------------------------------}

Function TransformadaAPlanoInclinado(ndia:dia;beta:real;modelo:string;Hg:real):real;
var
  hora: integer;
  global,directa,difusa,Htotal:real;
  wdif,wmed,cosho,cosinc: real;
  tmh: real;
  j: integer;
begin

  model:= -1;
		 if modelo='Liu-Jordan' then model:=1;
		 if modelo='HDKR' then model:=2;
		 if modelo='Perez' then model:=3;

{1}    Htotal:=0;
		 dias:=ndia;


		 if dias= 366 then mes:=12 else mes:= trunc(dias/30.45 +1);

     {2}    Busqueda(Hg,registro);

{3}    Salida_puesta(ti,tf,wii,wff,tii,tff);

			 for hora:= 0 to 23 do EnergiaHoraPtr^[hora]:=0;

{4}    for j:=1 to 14 do
       begin
        hora:= j+4;
        if (hora>=tii) and (hora<tff) then
        begin
         global:= Hg* registro.fraccionglobal[j] ;
         directa:= global*registro.fracciondirecta[j];
         difusa:=  global-directa;

			Whora(hora,wmed,wdif);
         Verificacion( ndia, hora,wmed,wdif,global,directa,difusa);

			cosho:=
			cos_l*cos_delta[dias]*cos(wmed)+ sin_l*sin_delta[dias];
         cosinc:=
			cos_lbeta*cos_delta[dias]*cos(wmed) + sin_lbeta*sin_delta[dias];

	 case model of
		3: tmh:= Perez(beta,wdif,cosinc,cosho,global,directa,difusa);
		2: tmh:= HDKR(cosinc,cosho,global,directa,difusa);
		1: tmh:=	Liu_Jordan(cosinc,cosho,global,directa,difusa);
   else
    raise Exception.Create(' Modelo no válido, IT4->TransformadaAPlanoInclinado ');
	 end;
  if tmh <  0  then
    raise Exception.Create( 'tmh < 0 , IT4->TransformadaAPlanoInclinado');

	 HTotal:= Htotal+ tmh;
	 EnergiaHoraPtr^[hora]:= tmh;
end;
end;

TransformadaAPlanoInclinado:= Htotal;
end;
{--------------------------------------------}

Procedure IT_Init( beta:real);
var
  i, j: integer;
  jgkt: integer;
begin

     if beta>pi/2 then beta:=pi/2;
     cos_l:= cos(latitud);
    sin_l:= sin(latitud);
    cos_beta:= cos(beta);
    sin_beta:= sin(beta);
	 sin_betac:= sqr(sin(beta/2))*sin(beta/2);
     cos_lbeta:= cos(beta+latitud);
    sin_lbeta:= sin(beta+latitud);
    for i:=1 to 366 do
     begin
			delta[i]:= pi*23.45/180*sin(2*pi*(284+i)/365.25);

      cos_delta[i]:=cos(delta[i]);
        sin_delta[i]:=sin(delta[i]);

        Isn[i]:=12/pi*Isc*(1+0.033*cos(2*pi*i/365.25));
      Ws[i]:= Acos(-Tan(latitud)*Tan(delta[i]));
      sin_ws[i]:= sin(ws[i]);

       Ho[i]:= 3.6/41.86* 2*Isn[i]*(
       cos_l*cos_delta[i]*sin_ws[i] +  Ws[i]*sin_l*sin_delta[i] );
       { [Ho]=  cal/cm2 }

     end;
 {Lectura de indices y asignacion}

  ref[1]:= 95;
  ref[2]:= ref[1]+43;
  ref[3]:= ref[2]+74;
  ref[4]:= ref[3]+86;
  ref[5]:= ref[4]+153;
  ref[6]:= ref[5]+113;
  ref[7]:= ref[6]+132;
  ref[8]:= ref[7]+105;
  ref[9]:= ref[8]+116;
  ref[10]:= ref[9]+122;
  ref[11]:= ref[10]+108;
	ref[12]:= ref[11]+133;
 for jgkt:=1 to 390 do  total[jgkt]:=0;
 mes:=1;

 seek(archivo,0);
 while not eof(archivo) do
 begin
    read(archivo,registro);
    j:= filepos(archivo);
     if j> ref[mes]  then  mes:=mes+1;
     kt:= registro.ktdato;
    grupokt:= trunc(kt/2/step+1) +30*(mes-1) ;
    total[grupokt]:= total[grupokt] +1;
    indicePtr^[grupokt,total[grupokt]]:= j ;
  end;

{Datos para el modelo de Perez et al.}

efe[1,1]:=  -0.196;  efe[1,2]:=   1.084;  efe[1,3]:=   -0.006;
efe[1,4]:=  -0.114;  efe[1,5]:=   0.180;  efe[1,6]:=   -0.019;
efe[2,1]:=   0.236;  efe[2,2]:=   0.519;  efe[2,3]:=   -0.180;
efe[2,4]:=  -0.011;  efe[2,5]:=   0.020;  efe[2,6]:=   -0.038;
efe[3,1]:=   0.454;  efe[3,2]:=   0.321;  efe[3,3]:=   -0.255;
efe[3,4]:=   0.072;  efe[3,5]:=  -0.098;  efe[3,6]:=   -0.046;
efe[4,1]:=   0.866;  efe[4,2]:=  -0.381;  efe[4,3]:=   -0.375;
efe[4,4]:=   0.203;  efe[4,5]:=  -0.403;  efe[4,6]:=   -0.049;
efe[5,1]:=   1.026;  efe[5,2]:=  -0.711;  efe[5,3]:=   -0.426;
efe[5,4]:=   0.273;  efe[5,5]:=  -0.602;  efe[5,6]:=   -0.061;
efe[6,1]:=   0.978;  efe[6,2]:=  -0.986;  efe[6,3]:=   -0.350;
efe[6,4]:=   0.280;  efe[6,5]:=  -0.915;  efe[6,6]:=   -0.024;
efe[7,1]:=   0.748;  efe[7,2]:=  -0.913;  efe[7,3]:=   -0.236;
efe[7,4]:=   0.173;  efe[7,5]:=  -1.045;  efe[7,6]:=    0.065;
efe[8,1]:=   0.318;  efe[8,2]:=  -0.757;  efe[8,3]:=    0.103;
efe[8,4]:=   0.062;  efe[8,5]:=  -1.698;  efe[8,6]:=    0.236;

end;

{--------------------------------------------}




function NuevoDia(
	ndia:dia;   { D¡a del a¤o. 1 = primero de enero }
	beta:real;       { Inclinaci¢n del panel [radianes] }
	modelo:string; {modelo de rad. solar, existen 3 modelos
						 'Liu-Jordan' 'Perez' y 'HDKR'}
	Hg: real
	 ):real; { [Wh/m2] }
var
	k: integer;
begin
	Hg:=Hg/1000/100/100; {[KWH/cm2]}
	Hg:= KWH_A_Cal(Hg); { [Cal/cm2] }
	Hg:= TransformadaAPlanoInclinado( ndia, beta, modelo, Hg );
	NuevoDia:= Cal_A_KWH(Hg)*1000*100*100;
	for k:= 0 to 23 do
		EnergiaHoraPtr^[k]:=Cal_A_KWH(EnergiaHoraPtr^[k])*1000*100*100;
end;



function NuevoDia_HgCalcm2(
	ndia:dia;   { D¡a del a¤o. 1 = primero de enero }
	beta:real;       { Inclinaci¢n del panel [radianes] }
	modelo:string; {modelo de rad. solar, existen 3 modelos
						 'Liu-Jordan' 'Perez' y 'HDKR'}
	Hg: real
	 ):real; { [Wh/m2] }
var
	k: integer;
begin
//	Hg:=Hg/1000/100/100; {[KWH/cm2]}
//	Hg:= KWH_A_Cal(Hg); { [Cal/cm2] }
	Hg:= TransformadaAPlanoInclinado( ndia, beta, modelo, Hg );
	result:= Cal_A_KWH(Hg)*1000*100*100;
	for k:= 0 to 23 do
		EnergiaHoraPtr^[k]:=Cal_A_KWH(EnergiaHoraPtr^[k])*1000*100*100;
end;

function GetExeDir: string;
var
  s: string;
begin
    s:= ParamStr( 0 );
    s:= ExtractFileDir( s );
    result:= s;
end;

procedure logwrite( s: string );
var
  f: textfile;
begin
  assignfile( f, 'log.txt' );
  {$I-}
  append( f );
  {$I+}
  if ioresult <> 0 then
    rewrite( f );
  writeln( f, s );
  closefile( f );
end;


procedure datos2xlt;
var
  r: registrodedatos;
  f: textfile;
  j: integer;
begin
  assignfile( f, 'adatos.xlt' );
  rewrite( f );
  reset( archivo );
  while not eof ( archivo ) do
  begin
    read( archivo, r );
    write( f, r.ktdato );
    for j:= 1 to 14 do
    begin
      write( f, #9, r.fraccionglobal[j] );
    end;
    for j:= 1 to 14 do
    begin
      write( f, #9, r.fracciondirecta[j] );
    end;
    writeln( f );
  end;
  closefile( f );
  reset( archivo );
end;
begin
  sysutils.DecimalSeparator:= '.';
  sysutils.DateSeparator:= '/';

	new(EnergiaHoraPtr);
	new(IndicePtr);

  filemode:= 0;
  assign(archivo, GetExeDir+'\adatos');
  {$I-}
  reset(archivo);
  {$I+}
  if ioresult <> 0 then
    raise Exception.Create('Error!!, uidad: it4,  IT_Init, no puede abrir archivo '''+GetExeDir+'\adatos''');

//  datos2xlt;
  
	IT_Init(45*pi/180);

  logwrite('hola-inicio');
end.


