//  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 [pts] = CL_getEarthMap(res, coord, er)
// Earth map outlines
//
// Calling Sequence
// pts = CL_getEarthMap([res, coord, er])
//
// Description
// <itemizedlist><listitem>
// <p>Returns Earth map outlines coordinates (the same maps as used by CL_plot_earthMap()).</p>
// <p>The result is a 3xN matrix (spherical or geodetic coordinates). 
// The third row is initialized to 0 or "er" depending on the type of coordinates 
// ("er" for "sph", 0 for "ell"). </p>
// <p>The data consists in continuous segments separated by %nan. </p>
// <p></p>
// <p></p></listitem>
// <listitem>
// <p>Notes: </p>
// <p>- Earth map files are loaded from the directory: CL_home()/data/earthMap. </p>
// <p>- The file format is specific to CelestLab.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// res: (string, optional) Resolution: "low" or "high". Default is "low". 
// coord: (string, optional) Type of coordinates: "sph"=spherical, "ell"=elliptical. Default is "sph". 
// er: (optional) Equatorial radius [m]. Default value is [] (internal value used).
// pts: [longitude; latitude; z]; z = "er" if coord="sph", z = 0 if coord="ell". (3xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Loads a map (spherical coordinates, default radius)
// pts = CL_getEarthMap(); 
// // Plots the map
// CL_plot_ephem(CL_co_sph2car(pts)); 


// Notes: 
// The map is stored as a matrix of 16bits integers to save disk space. 
// Accuracy is about 6.e-3 degrees
// (%nan between blocks)
// 
// NB: 
// %CL_eqRad is not used as this function is for Earth only (%CL_eqRad is generic). 


if (~exists("coord","local")); coord = "sph"; end 
if (~exists("res","local")); res = "low"; end 
if (~exists("er","local")); er =[]; end 

if (er == []); er = CL__dataGetEnv("eqRad"); end


// file path
fname = fullfile(CL_home(), "data", "earthMap", "earthMap_" + coord + "_" + res + ".dat");

if (~isfile(fname))
  CL__error("File: " + fname + ": not found"); 
end

load(fname, "M"); 

// %nan values are stored in the 16bits integer matrix as -2^15
NAN = int16(-(2^15)); 
Inan = find(M(1,:) == NAN | M(2,:) == NAN);
  
// convert to double (degrees)
F = 2^15 - 1;
M = double(M) * 180 / F;

// z = value of 3rd coordinate
if (coord == "sph")
  z = er * ones(M(1,:)); 
else
  z = zeros(M(1,:)); 
end

// final results - radians (with %nan values)
M(1:2,Inan) = %nan;
z(Inan) = %nan;

pts = [ M * %pi/180; z ]; 
      
endfunction