{+doc
+NOMBRE: autoesca
+CREACION: 1.1.1990
+AUTORES: rch
+REGISTRO:
+TIPO: Unidad Pascal.
+PROPOSITO:Sevicios de calculo de escalas automaticas
+PROYECTO:rchlib

+REVISION: mayo 2017 rch, le agrego EscalaFechaN para graficos con eje
  fecha.

Busca la mejor aproximación del rango a un múltiplo
de  Año, dia, h, min, s

+AUTOR:
+DESCRIPCION:
-doc}

unit AutoEsca;

interface

uses
  Math, xMatDefs;

type

  TEscala = class
    grid, n, a: integer;
    xm1, xm2:   NReal;
    constructor Create(x1, x2: NReal);
  end;

  // dado xx1, xx2 y N calcula la escala que mejor aproxima con intervalos 10^b
  // siendo b = 1, 2 o 5. Si MODO = 0 aproxima por dentro si MODO = 1 aproxima por fuera.
procedure Escala125( out x1, x2, dx: NReal; xx1, xx2: NReal;  N: integer; MODO: integer );

procedure Escala125N( var x1, x2, dx: NReal; var N: integer; MODO: integer);
{
  Entradas:
    x1,x2 = rango  de valores que se desea cubrir.
    N = número de divisiones aprox. que se desea.
    MODO = modo en que se aproximará  el intervalo.
     MODO = 0, aproximar por adentro.
     MODO = 1, aproximar por afuera.
  Salidas:
    x1,x2 = rango que se cubrir . Calculado en la grilla de salida.
    Si Modo:=1 si el rango de salida aproxima al de entrada
    incluyendolo, por el contrario si Modo=0 lo aproxima sin
    incluir los extremos originales.(Lo aproxima por adentro).
    dx = ancho de una divisi¢n. Se cumple que:
      dx = an * exp10(a), donde an =(1,2,5) y a:integer.
    N = n£mero de divisiones que se tendr n.
}

procedure EscalaFechaN(var x1, x2, dx: NReal; var N: integer; MODO: integer);


implementation




function IntInf(x: NReal): int64;
begin
 result:= floor( x );
end;

function IntSup(x: NReal): int64;
begin
  result:= ceil( x );
end;


function Log(x: NReal): NReal;
begin
  result := log10( x );
end;

function Exp10(x: NReal): NReal;
begin
  result :=power( 10, x );
end;


// dado un rango x1 < x2, un número b > 0 y un número N de divisiones dados calcula el mínimo
// valor de k (entero) tal que existe un entero m tal que
//  ( m b 10^k <= x1 ) y ( x2<= (m+N) b 10^k  )
// la función retorna ( N b 10^k ) como resultado que es el ancho cubierto y
// los valores k y m por variable.
// Probando con b = 1, 2 y 5 se logra encontrar la mejor escala que aproxima
// por fuera un intervalo con N divisioines.
function min_Nb10k( var m, k: integer; x1, x2: NReal; b, N: integer): NReal;
var
  buscando: boolean;
  b10k: NReal;
  y2: NReal;

begin
  k:= IntSup( Log( (x2-x1) / ( N*b ) ));
  buscando:= true;
  while buscando do
  begin
    b10k:= b * Exp10( k );
    m:= intInf( x1 / b10k );
    y2:= ( m + N ) * b10k;
    if y2 >= x2 then
      buscando:= false
    else
      inc( k );
  end;
  result:= N * b10k;
end;



// dado un rango x1 < x2, un número b > 0 y un número N de divisiones dados calcula el mánimo
// valor de k (entero) tal que existe un entero m tal que
//  ( x1 <= m b 10^k ) y ( (m+N) b 10^k  <= x2  )
// la función retorna ( N b 10^k ) como resultado que es el ancho cubierto y
// los valores k y m por variable.
// Probando con b = 1, 2 y 5 se logra encontrar la mejor escala que aproxima
// por dentro un intervalo con N divisioines.
function max_Nb10k( out m, k: integer; x1, x2: NReal; b, N: integer): NReal;
var
  buscando: boolean;
  b10k: NReal;
  y2: NReal;

begin
  k:= IntSup( Log( (x2-x1) / ( N*b ) ));
  buscando:= true;
  while buscando do
  begin
    b10k:= b * Exp10( k );
    m:= intSup( x1 / b10k );
    y2:= ( m + N ) * b10k;
    if y2 <= x2 then
      buscando:= false
    else
      dec( k );
  end;
  result:= N * b10k;
end;




constructor TEscala.Create(x1, x2: NReal);
var
  dx:   NReal;
  temp: integer;

  procedure v(gridx, nx: integer);
  begin
    n    := nx;
    grid := gridx;
  end;

begin
  inherited Create;

  dx   := x2 - x1;
  a    := IntSup(log(dx / 50));
  temp := IntSup(dx * exp10(-a));
  case temp of
    5: v(5, 1);
    6..8: v(8, 1);
    9..10: v(10, 1);
    11..16: v(8, 2);
    17..20: v(10, 2);
    21..25: v(5, 5);
    26..40: v(8, 5);
    41..50: v(10, 5);
  end;
  temp := IntInf(x1 / exp10(a) / n);
  xm1  := temp * exp10(a) * n;
  xm2  := (temp + grid) * exp10(a) * n;
end;



procedure Escala125(out x1, x2, dx: NReal; xx1, xx2: NReal; N: integer;
  MODO: integer);
var
  b: array[1..4] of integer;
  m, k: Integer;
  Nb10k, Nb10k_opt: NReal;
  j: integer;
  aux: NReal;
  tolerancia_dx, xx1_, xx2_: NReal;
begin
    if xx1 > xx2 then
    begin // hay que ayudar a los tontos
      xx1_:= xx2;
      xx2_:= xx1;
    end
    else
    begin
      xx1_:= xx1;
      xx2_:= xx2;
    end;

    // altermaos un poquito los extremos para que no genere una nueva división
    // si estamos en la tolerancia
    tolerancia_dx:= ( xx2_ - xx1_ ) / N / 1000;
    if tolerancia_dx < 1e-10 then
      tolerancia_dx:= 1e-10;

    if MODO = 0 then
    begin
      xx1_:= xx1_ - tolerancia_dx;
      xx2_:= xx2_ + tolerancia_dx;
    end
    else
    begin
      xx1_:= xx1_ + tolerancia_dx;
      xx2_:= xx2_ - tolerancia_dx;
      if xx1_ > xx2_ then
      begin
        aux:= xx1_;
        xx1_:= xx2;
        xx2_:= aux;
      end;
    end;

    if ( xx2_ - xx1_ ) < (2*tolerancia_dx) then
       xx2_:= xx1_ + 2* tolerancia_dx;


  b[1]:= 1;
  b[2]:= 2;
  b[3]:= 5;
  b[4]:= 25;

  if MODO = 0 then
  begin
    Nb10k_opt:= 0; // seguro esto es cota infior de los 3
    for j:= 1 to 4 do
    begin
      Nb10k:= max_Nb10k( m, k, xx1_, xx2_, b[j], N );
      if Nb10k > Nb10k_opt then
      begin
        dx:= Nb10k/N;
        x1:= m * dx;
        x2:= x1 + Nb10k;
        Nb10k_opt:= Nb10k;
      end;
    end;
  end
  else
  begin
    Nb10k_opt:= 100* ( xx2 - xx1 )+100;  // seguro esto es cota superior
    for j:= 1 to 4 do
    begin
      Nb10k:= min_Nb10k( m, k, xx1_, xx2_, b[j], N );
      if Nb10k < Nb10k_opt then
      begin
        dx:= Nb10k/N;
        x1:= m * dx;
        x2:= x1 + Nb10k;
        Nb10k_opt:= Nb10k;
      end;
    end;
  end
end;

procedure Escala125N(var x1, x2, dx: NReal; var N: integer; MODO: integer);
var
  y1, y2, dy: array[-2..2] of NReal;
  j: integer;
  dy1y2, dy1y2_opt: NReal;
  j_opt: integer;
begin

  if Modo = 1 then
  begin
    dy1y2_opt := (x2-x1)*100+100;
    for j := -2 to 2 do
    begin
      if ( N+j )> 0 then
      begin
        Escala125( y1[j], y2[j], dy[j], x1, x2, N+j, 1 );
        dy1y2:= ( y2[j]- y1[j] );
        if dy1y2  < dy1y2_opt then
        begin
          j_opt:= j;
          dy1y2_opt:= ( y2[j]- y1[j] );
        end;
      end;
    end;
  end
  else
  begin
    dy1y2_opt := 0;
    for j := -2 to 2 do
    begin
      if ( N+j )> 0 then
      begin
        Escala125( y1[j], y2[j], dy[j], x1, x2, N+j, 0 );
        dy1y2:= ( y2[j]- y1[j] );
        if dy1y2  > dy1y2_opt then
        begin
          j_opt:= j;
          dy1y2_opt:= ( y2[j]- y1[j] );
        end;
      end;
    end;
  end;
  x1:= y1[j_opt];
  x2:= y2[j_opt];
  dx:= dy[j_opt];
  N:= N + j_opt;

end;

     (*
procedure Escala125N(var x1, x2, dx: NReal; var N: integer; MODO: integer);
var
  m, aux: NReal;
  n1, n2: int64;
  k, j: integer;
  temp:   array[1..3] of NReal;
  a:      array[1..3] of integer;
  expo:   NReal;
begin
  if x1 > x2 then
  begin // hay que ayudar a los tontos
    aux:= x1;
    x1:= x2;
    x2:= aux;
  end;

  if ( x2 - x1 ) < 1E-32 then
  begin
    if abs( x1 ) > 1E-10 then
      x2:= x1 * 1.1
    else
      x2:= 1E-10;
    if x1 > x2 then
    begin
      aux:= x1;
      x1:= x2;
      x2:= aux;
    end;
  end;

  aux:= log( x2 - x1 );

  temp[1] := log( N );
  temp[2] := log(2 * N);
  temp[3] := log(5 * N);

  if MODO = 1 then
  begin
    m := 1e30;
    j := 0;
    for k := 1 to 3 do
    begin
      a[k] := IntSup(aux - temp[k]);
      if temp[k] + a[k] < m then
      begin
        m := temp[k] + a[k];
        j := k;
      end;
    end;
  end
  else
  begin
    m := -1e30;
    j := 0;
    for k := 1 to 3 do
    begin
      a[k] := IntInf(aux - temp[k]);
      if temp[k] + a[k] > m then
      begin
        m := temp[k] + a[k];
        j := k;
      end;
    end;
  end;

  if j = 3 then
    dx := 5
  else
    dx := j;

  expo := a[j];
  aux  := exp10(expo) * dx;
  if MODO = 1 then
  begin
    n1 := IntInf{Sup}(x1 / aux);
    n2 := IntSup{Inf}(x2 / aux);
  end
  else
  begin
    n1 := IntSup(x1 / aux);
    n2 := IntInf(x2 / aux);
  end;
  n  := n2 - n1;
  x1 := n1 * aux;
  x2 := n2 * aux;
  dx := aux;
end;

*)



function EnRango20( x: NReal ): boolean;
begin
  result:= ( x > 0.8 ) and ( x < 1.2 );
end;



const
  dtAnio = (97.0 * 366.0 + (400.0 - 97.0) * 365) / 400.0;
  dtMes = dtAnio  / 12;
  dtTrimestre = dtAnio / 4;
  dtSemestre  = dtAnio / 2;

  divisores_dia: array[0..12] of NReal = (
1.0/24.0/60.0 , // 0
1.0/24.0 ,      // 1
3.0/24.0 ,      // 2
4.0/24.0 ,      // 3
6.0/24.0 ,      // 4
12.0/24.0 ,      // 5
1.0,            // 6
7.0,            // 7
15.0,           // 8
dtMes,           // 9
dtTrimestre,           // 10
dtSemestre,      //11
dtAnio );        // 12


procedure EscalaFechaN(var x1, x2, dx: NReal; var N: integer; MODO: integer);
var
  nDivisores: integer;
  a_q, best_q: NReal;
  d1, d1_best: NReal;
  k_best_q: integer;
  k: integer;

begin
  if x2 < x1 then vswap( x1, x2 );
  dx:= ( x2-x1 ) / N;

  nDivisores:= length( divisores_dia );
  k_best_q:= 0;
  best_q:= dx/divisores_dia[0];
  d1_best:= abs( 1 - best_q );
  for k:= 1 to NDivisores - 1 do
  begin
    a_q:= dx/divisores_dia[k];
    d1:= abs( 1 - a_q );
    if d1 < d1_best then
    begin
      k_best_q:= k;
      d1_best:= d1;
      best_q:= a_q;
    end;
  end;


  dx:= divisores_dia[ k_best_q ];
  if best_q > 1 then
    dx:= dx * round( best_q );
  if k_best_q >= 6 then
    dx:= round( dx );

  if modo = 0 then
  begin // por dentro
    x1:= IntSup( x1 / dx )* dx;
    x2:= IntInf( x2 / dx )* dx;
  end
  else
  begin // por fuera
    x1:= IntInf( x1 / dx )* dx;
    x2:= IntSup( x2 / dx )* dx;
  end;
  N:= Round( ( x2 - x1 ) / dx );

end;


end.

