{+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;

interface

type
  dia = 1..366;


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

var
  EnergiaHoraPtr: PVectDatHr;


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

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

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

function cal_A_kWh(cal: double): double;
function kWh_A_cal(KWH: double): double;

procedure IT_Init(beta: double);


implementation


const
  pi = 3.141592;   {n£mero pi}
  l = -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 double;

  vector2 = array[dia] of double;


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

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

  Isn, delta, Ws, Ho, cos_Ws, 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 double;

  tii, tff, j, grupokt: integer;

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

  model: integer;

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

function CALporCM2_A_KJporM2(x: double): double;
begin
  CALporCM2_A_KJporM2 := x * 41.86;
end;

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

function Tan(x: double): double;
begin
  Tan := sin(x) / cos(x);
end;

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

function Acos(x: double): double;
begin
  Acos := -Arctan(x / sqrt(1 - x * x)) + pi / 2;
end;

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

function cal_A_kWh(cal: double): double;
begin
  cal_A_kWh := cal * 1.162E-6;
end;

function kWh_A_cal(kWh: double): double;
begin
  kWh_A_cal := kWh / 1.162E-6;
end;

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

function Ktf(Hg: double): double;
begin
  Ktf := Hg / Ho[dias];
end;

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

procedure Busqueda(Hg: double; 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: double; var tii, tff: integer);
var
  td: double;
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: double);

var
  w1, w2: double;

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(hora, wmed, wdif: double; var global, directa, difusa: double);
var
  kthora: double;

begin
  I0hora := 3.6 / 41.86 * Isn[i] * (cos_l * cos_delta[i] * 2 *
    sin(wdif / 2) * cos(wmed) + wdif * sin_l * sin_delta[i]); { cal/cm2}
  kthora := global / I0hora;
  if kthora > 1 then
  begin
    writeln('corrijo global dias hora  ', dias, '   ', hora);
    directa := directa / global * 0.9 * I0hora;
    global := 0.9 * I0hora;
    difusa := global - difusa;
  end;
end;

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

function Liu_Jordan(cosinc, cosho, global, directa, difusa: double): double;

var
  reflejada: double;
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: double): double;
  {HAY DAVIES KLUCHER REINDL}

var
  f, ai, reflejada: double;

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: double): double;

var

  aa, bb, ab, exn, del, directan, titaz, eps, f1, f2, tit3: double;
  reflejada: double;
  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 <> 0 then
    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];

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

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

  Perez := directa + difusa + reflejada;
end;

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

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

begin

  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(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);
      end;
      HTotal := Htotal + tmh;
      EnergiaHoraPtr^[hora] := tmh;
    end;
  end;

  TransformadaAPlanoInclinado := Htotal;
end;

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

procedure IT_Init(beta: double);
begin

  writeln('Randomize');

  Randomize;

  writeln('VolvĦ del randomize');

  if beta > pi / 2 then
    beta := pi / 2;
  cos_l := cos(l);
  sin_l := sin(l);
  cos_beta := cos(beta);
  sin_beta := sin(beta);
  sin_betac := sqr(sin(beta / 2)) * sin(beta / 2);
  cos_lbeta := cos(beta + l);
  sin_lbeta := sin(beta + l);
  writeln('hola1');
  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(L) * 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;
  writeln('hola2');
  {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;
  writeln('hola3');
  for grupokt := 1 to 390 do
    total[grupokt] := 0;
  mes := 1;

  writeln('hola4');
  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;

  writeln('hola5');
  {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: double;       { Inclinaci˘n del panel [radianes] }
  modelo: string; {modelo de rad. solar, existen 3 modelos
             'Liu-Jordan' 'Perez' y 'HDKR'}
  Hg: double): double; { [cal/cm2] }
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;




begin
  writeln('it4-begin');
  writeln('new(EnergiaHoraPtr)');
  new(EnergiaHoraPtr);
  writeln('new(IndicePtr)');
  new(IndicePtr);

  Assign(archivo, 'adatos.r8');
        {$I-}
  reset(archivo);
        {$I+}
  if ioresult <> 0 then
  begin
    writeln('Error!!, uidad: it4,  IT_Init, no puede abrir archivo ''adatos''');
    writeln('presione ENTER para finalizar..');
    readln;
    halt(1);
  end;


  IT_Init(45 * pi / 180);
  writeln('it4-end');
end.
