//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.
 
function [val] = CL_locationInfo(info,lon,lat)
// Information about a location on the ground (land ratio)
//
// Calling Sequence
// val = CL_locationInfo("land_ratio",lon,lat)
//
// Description
// <itemizedlist><listitem>
// <p>Information about a location on the ground.</p>
// <p></p>
// <p>Only option available: "land_ratio", which gives the percentage 
// of land in a 0.25 deg^2 area around the location: </p>
// <p>- 1.0 : 100% land, 0% sea </p>
// <p>- 0.7 : 70% land, 30% sea </p>
// <p>- 0.0 : 0% land, 100% sea </p>
// </listitem>
//</itemizedlist>
//
// Parameters
// info: (string) Type of information to be returned (1x1)
// lon : Geodetic longitude of the location [rad] (PxN)
// lat: Geodetic latitude of the location [rad] (PxN)
// val : Information value returned for each location (PxN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// lon = CL_deg2rad(-10: 0.5: 40);
// lat = CL_deg2rad(30: 0.5: 60); 
// [lon1,lat1] = ndgrid(lon,lat); 
// land_ratio = CL_locationInfo("land_ratio", lon1, lat1);
//
// // Plot results on an Earth map:
// f = scf();
// f.color_map = jetcolormap(50);
// Sgrayplot(CL_rad2deg(lon), CL_rad2deg(lat), land_ratio, colminmax=[1,50]);
// CL_plot_earthMap(color_id = color("black"), data_bounds=[-10,30;40,60], ..
// res = "high", coord="ell", tick_steps=[10,10]);


function [desc,data] = read_info(id)
// Loads the info from a binary file
// NB: the file contains the values for each grid point
// The set of points is described by "desc" (nearly uniformly 
// distributed points in surface). 
// The file contains the size of the mesh (number of latitudes). 
// and the data for all the mesh (as a matrix of bytes to save memory)
// id: (string) info identifier
// desc: structure describing the meshing
// data: (1xN) matrix of doubles: value in [0,1] 

  // file path
  fname = fullfile(CL_home(), "data", "earthMap", "info_" + id + ".dat");
  
  if (~isfile(fname))
    CL__error(id + " not recognized as a valid type of information");
  end
  
  // Load the binary file that contains variables "nblat" and "info"
  load(fname, "nblat", "info");
  
  // Creation of structure describing mesh
  desc = gr2d_new(nblat);
  
  // Convert info from short integers : [0-255] to double [0-1]
  F = 2^8 - 1; 
  data = double(info) / F;
     
endfunction

// -------------------------------------------------
// Creation of structure describing mesh
// -------------------------------------------------
function [desc] = gr2d_new(nblat)
  if (nblat < 1 | round(nblat) <> nblat); CL__error("*** Invalid number"); end

  paslat = %pi / nblat; 
  latmin = -%pi/2 + paslat/2; 
  latmax =  %pi/2 - paslat/2; 
  tab_lat = linspace(latmin, latmax, nblat); 
  nblon = round( 2*%pi * cos(tab_lat) ./ paslat ) ;
  paslon = 2*%pi ./ nblon; 
  lonmin = -%pi + paslon/2; 

  desc = struct();
  desc.ynb = nblat; 
  desc.ymin = latmin; 
  desc.yp = paslat; 
  desc.xnb = nblon; 
  desc.xmin = lonmin;
  desc.xp = paslon; 
endfunction


// -------------------------------------------------
// Interpolation:
// -------------------------------------------------
function [z1] = gr2d_interp(desc, z, x1, y1)
  dims1 = size(x1);
  x1 = matrix(x1(:),1,-1); // ligne 
  y1 = matrix(y1(:),1,-1); // ligne
  
  z = matrix(z, 1, -1); 

  D = desc;
  D.xind = [1, 1+cumsum(D.xnb(1:$-1))]; 
  
  deff("[xij] = X(i,j)", "xij = D.xmin(j)+(i-1).*D.xp(j)"); 
  deff("[yj] = Y(j)", "yj = D.ymin+(j-1)*D.yp"); 
  deff("[zij] = Z(i,j)", "zij = z(D.xind(j) + i - 1)"); 

  j = min(max(floor((y1-D.ymin)/D.yp)+1,1),D.ynb-1); 
  i1 = min(max(floor((x1-D.xmin(j))./D.xp(j))+1,1),D.xnb(j)-1); 
  i2 = min(max(floor((x1-D.xmin(j+1))./D.xp(j+1))+1,1),D.xnb(j+1)-1); 

  alpha1 = ((X(i2, j+1) - X(i1,j)) ./ D.xp(j)); 
  alpha2 = ((X(i2+1, j+1) - X(i1,j)) ./ D.xp(j)); 

  a = Z(i1,j); 
  b = Z(i1+1,j) - a;  
  d = (Z(i2+1,j+1) - Z(i2,j+1)) ./ (alpha2 - alpha1) - b; 
  c = Z(i2,j+1) - a - (b+d) .* alpha1; 

  x2 = ((x1 - X(i1,j)) ./ D.xp(j)); 
  y2 = ((y1 - Y(j)) / D.yp); 

  z1 = matrix(a + b .* x2 + c .* y2 + d .* x2 .* y2, dims1); 
endfunction


// -------------------------------------------------
// Main
// -------------------------------------------------

// Declarations:


// Code: 

// error checking
if (size(lon,1) ~=  size(lat,1) | size(lon,2) ~=  size(lat,2))
  CL__error("lon and lat must have the same size"); 
end


// loading of the info:
[desc,data] = read_info(info);

// set longitude between -%pi and %pi
lon = CL_rMod(lon,-%pi,%pi);

// interpolate
val = gr2d_interp(desc, data, lon, lat);

endfunction