//  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'.

// -----------------------------------------------------------
//> Earth's magnetic field according to IGRF-11 model. 
//> The parameters that can be plotted are: 
//> -> Components of magnetic field: 
//> North-South (+=N), East-West (+=E), Vertical (+=descending) 
//> -> Intensity of magnetic field (>=0): 
//> horizontal or total 
//> -> Direction of magnetic field: 
//> inclination / horizontal plane (+=descending)
//> declination / North (+=East) 
//
// Author: A. Lamy
// -----------------------------------------------------------

// ---------------------------------------------------------
// Check if dates are in the correct range
// NB: cal is invalid if it contains %nan (then => %f)
// scal: calendar string, dur = [d_past, d_future] (days)
// ---------------------------------------------------------
function [ok] = isValidDate(scal)  
  ok = %f; 
  cjd1 = CL_dat_cal2cjd(1900,1,1);
  cjd2 = CL_dat_cal2cjd(2020,1,1);
  cal = CL_dat_str2cal(scal); 
  if (cal(1) <> %nan)
    cjd = CL_dat_cal2cjd(cal); 
    ok = ( cjd >= cjd1 & cjd <= cjd2 );
  end
endfunction


// -----------------------------------------------------------
// Parameters input
// -----------------------------------------------------------

// Initial start date = current time (TREF)
scal = CL_dat_cal2str(CL_dat_cjd2cal(floor(CL_dat_now()))); 
alt = 0; // geodetic altitude 
opt = 4; 

desc = list(..
  CL_defParam("Date (calendar format, TREF time scale)", scal, typ='cal', valid="isValidDate($x)"), ...
  CL_defParam("Geodetic altitude", alt, units=["m", "km"], valid="$x >= 0 & $x <= 50000"), ...
  CL_defParam("Option (1=comp-N, 2=comp-E, 3=comp-V, 4=int_hor, 5=int_tot, 6=incl, 7=decl)", opt, accv=1:7) ...
);

[scal, alt, opt] = CL_inputParam(desc);


// -----------------------------------------------------------
// computations
// -----------------------------------------------------------

cjd = CL_dat_cal2cjd(CL_dat_str2cal(scal));

// position - geodetic coordinates; 
lon = linspace(-%pi, %pi, 80); 
lat = linspace(-%pi/2, %pi/2, 60); 

[Lon, Lat] = ndgrid(lon,lat); 

pos_ell = [matrix(Lon,1,-1); matrix(Lat,1,-1); alt * ones(matrix(Lon,1,-1))]; 

pos = CL_co_ell2car(pos_ell); 

// magnetic field
magf = CL_mod_geomagField(cjd, pos); 

// field in local frame (topocentric North)
magf_topo = CL_fr_topoNMat(pos_ell) * magf; 

// components: north, east, descending
// intensity
// direction
magf_n   = matrix(magf_topo(1,:), size(Lon)); 
magf_e   = matrix(-magf_topo(2,:), size(Lon)); 
magf_d   = matrix(-magf_topo(3,:), size(Lon)); 
magf_hor = matrix(CL_norm(magf_topo(1:2,:)), size(Lon)); 
magf_tot = matrix(CL_norm(magf_topo(1:3,:)), size(Lon)); 
magf_inc = atan(magf_d, magf_hor); 
magf_dec = atan(magf_e, magf_n); 


// -----------------------------------------------------------
// plot
// -----------------------------------------------------------

f = scf(); 
f.visible="on";
f.immediate_drawing="off";

nbcol = 128;
cmap = rainbowcolormap(nbcol); 
f.color_map = 0.2 + 0.8 * cmap($:-1:1,:);


// parameter to be plotted (nanotesla or deg)
select (opt)
  case 1; z = magf_n * 1.e9; 
  case 2; z = magf_e * 1.e9; 
  case 3; z = magf_d * 1.e9; 
  case 4; z = magf_hor * 1.e9; 
  case 5; z = magf_tot * 1.e9; 
  case 6; z = magf_inc * %CL_rad2deg; 
  else;   z = magf_dec * %CL_rad2deg; 
end

zmin = min(z); 
zmax = max(z); 

colorbar(zmin,zmax,colminmax=[1,nbcol],fmt="%.1f");

Sgrayplot(lon*%CL_rad2deg, lat*%CL_rad2deg, z, zminmax=[zmin, zmax], colminmax=[1,nbcol]);

// levels
zmin1 = CL_quantile(z(:)', 0.1); 
zmax1 = CL_quantile(z(:)', 0.9); 

levels = CL_autoLevels(zmin1, zmax1, 4);
// add levels
delta = levels($) - levels(1); 
levels = [levels-delta, levels, levels	+delta];
contour2d(lon*%CL_rad2deg, lat*%CL_rad2deg, z, levels, style=ones(levels)*color("grey30"));
h = CL_g_select(gce(),"Polyline"); 
h.thickness = 1;
h = CL_g_select(gce(), "Text");
h.font_foreground = color("grey10");
h.font_size = 1;
CL_g_set(h, "text", string(strtod(h.text)));


CL_plot_earthMap(color_id=color("grey10"), coord="ell"); 


// titles

a=gca();

txt = [["North component", "East component", "vert. component (+=D)", "horiz. intensity", "total intensity"] + " (nanotesla)", ["inclination (+=D)", "declination (+=E)"] + " (deg)"]; 
a.title.text = "Magnetic field " + txt(opt) + " - Alt = " + sprintf("%.0f", alt/1000) + " km"; 
a.x_label.text = "Geodetic longitude"; 
a.y_label.text = "Geodetic latitude"; 


f.immediate_drawing="on";
f.visible="on";