function PlotMapBougainville(img_name, subset_num, maptype, variable_toplot)
% plot satellite data with surface currents, with EEZ border and 12nm limits, and removing shallow pixels
% Author: Guillaume Bourdin
% Date: 2023-10-16
%% define path of folder containing satellite image and mercator netCDF file
addpath(genpath('/Users/gui/Documents/Maine/Data/TaraPacific/remote_sensing/time_series/level3_timeseries/Bougainville_MapPackage'))
addpath(genpath('/Users/gui/Documents/Maine/Data/TaraPacific/remote_sensing/time_series/level3_timeseries/highres_IMEdetection'))
path_bathymetry = '/Users/gui/Documents/Maine/Data/TaraPacific/remote_sensing/time_series/level3_timeseries/highres_IMEdetection/data/GEBCO_2022_sub_ice_topo.nc';
isl_id = img_name(1:3);
path_containg_img = ['/Users/gui/Documents/Maine/Data/TaraPacific/remote_sensing/time_series/level3_timeseries/' isl_id];

% set email details
emailfrom = 'satellite.bougainville@gmail.com';
pwd = 'zmzm nkfv dehc gyyp';
emailto = {'guillaume.bourdin@maine.edu'};
% emailto = {};
% emailto = {'vargas@sb-roscoff.fr', 'emmanuel.boss@maine.edu', 'fabien.lombard@imev-mer.fr', ...
%   'noan@planktonplanet.org', 'nils.haentjens@maine.edu', 'christophe.prazuck@sorbonne-universite.fr',...
%   'annaodd@gmail.com'};

if ~strcmp(img_name(end-3:end), '.mat')
  img_name = [img_name '.mat'];
end
% load satellite level-3 composite 
load(fullfile(path_containg_img, img_name), 'lat', 'lon', 'data_out', 'selOC')

mld_bool = false;
% Make sure all data come in the right shape: latitude gradient over dimension 1 and longitude on dimension 2
if abs(mean(diff(lat(1,:)))) > abs(mean(diff(lat(:,1))))
  fprintf('Latitude gradient detected on dim 2 and lon on dimension 1, transposing %s data ... ', fullfile(path_containg_img, img_name))
  lat = lat';
  lon = lon';
  for tr = 1:size(data_out,1)
    data_out.data_avg{tr} = data_out.data_avg{tr}';
    data_out.data_std{tr} = data_out.data_std{tr}';
    data_out.data_n{tr} = data_out.data_n{tr}';
  end
end
%% %%%%%%%%%%%%%%% Initiate station and email info %%%%%%%%%%%%%%% %%
if any(strcmp(isl_id, {'brr', 'brt', 'brm', 'bre'}))
  station = table();
  station.name = {'    Walters shoal'; '    MAD-Ridge'; '    La Réunion';
    '    Mauritius'; '    Tromelin'; '    Les Glorieuses'; '    Dar es Salaam'; 
    '    Mayotte'; '    Juan de Nova'; '    Bassas da India'; '    Europa'};
  station.lat = [-33.20; -27.5; -21.1514; -20.2863; -15.8913; -11.5801; -6.8246; -12.7882; -17.0564; -21.4754; -22.3805];
  station.lon = [43.8333; 46.2667; 55.5455; 57.5835; 54.5230; 47.2957; 39.2848; 45.2375; 42.7239; 39.6949; 40.3580];
  % emailto = [emailto {'manon.thueux@intradef.gouv.fr', 'thomas-f.finet@intradef.gouv.fr'}];
  emailto = [emailto {'guillaume.bourdin@maine.edu'}];
  voab = 'Manon et Thomas';
elseif any(strcmp(isl_id, {'bcv', 'bcs', 'bcp', 'bcc', 'bcf'}))
  station = table();
  station.name = {'    Nouméa'; '    Port Villa (Vanuatu)'; '    Vanuatu'; ...
    '    Nendo'; '    Solomon Isl.'; '    Fiji'; '    Tonga';};
  station.lat = [-22.247496; -17.734083; -15.4647; -10.7146; -9.809726; -17.828937; -20.247718];
  station.lon = [166.457909; 168.313611; 166.8497; 165.9359; 160.638263; 178.038627; -174.759111];
  % emailto = [emailto {'mathilde-a.vigneron@intradef.gouv.fr', 'hugo.zaccomer@intradef.gouv.fr'}];
  emailto = [emailto {'guillaume.bourdin@maine.edu'}];
  voab = 'Mathilde et Hugo';
elseif any(strcmp(isl_id, {'bpt', 'bps', 'bpm', 'bpg', 'bpc'}))
  station = table();
  station.name = {'    Papeete'; '    Raiatea'; '    Rurutu'; ...
    '    Gambier'; '    Rarotonga (Cook)'; '    Samoa'; '    Fakarava (Tuamotu)'; ...
    '    Nuku Hiva (Marquises)'; '    Clipperton'; '    Tonga';};
  station.lat = [-17.539551; -16.769091; -22.473030; -23.117973; -21.235415; -13.787905; -16.095449; -8.867158; 10.303843; -20.247718]; 
  station.lon = [-149.580404; -151.455468; -151.336651; -134.972668; -159.777030; -172.117330; -145.654007; -140.120811; -109.217633; -174.759111];
  % emailto = [emailto {'mathilde-a.vigneron@intradef.gouv.fr', 'hugo.zaccomer@intradef.gouv.fr'}];
  emailto = [emailto {'guillaume.bourdin@maine.edu'}];
  voab = 'Emmanuel et Nils';
else
  error('Leg ID unknown')
end

%% wrap longitude around 360°
lon360 = lon;
if all(any(lon(:) < 0) & any(lon(:) >= 0) & ...
    sum(lon(:) < -90 | lon(:) > 90) > sum(lon(:) > -90 & lon(:) < 90))
  lon360(lon360 < 0) = lon360(lon360 < 0) + 360;
end
full_img = [min(lat(:)) min(lon360(:)) max(lat(:)) max(lon360(:))];
if full_img(2) > 180
  full_img(2) = full_img(2) - 360;
end
if full_img(4) > 180
  full_img(4) = full_img(4) - 360;
end

%% %%%%%%%%%%%%%%% Initiate image subest depending on leg ID %%%%%%%%%%%%%%% %%
subset = table();
email_text = table();
switch isl_id
  case 'brr'
    subset.name = {'img_entiere'; 'Reunion_Mauritius'; 'Tromelin';'Sud_Reunion'};
    subset.limits = [full_img; % #1 img_entiere
      -15.2627 58.9733 -23.4285 50.2862; % #2 Reunion/Mauritius
      -17.6914 52.4069 -13.5306 56.5576; % #3 Tromelin
      -28.0262 46.1326 -19.8447 58.5919]; % #4 Reunion Sud
    email_text.sbj = {'Cartes Bougainville - La Reunion';
      'zoom La Reunion - Mauritius';
      'zoom Tromelin';
      'Sud La Reunion'};
    email_text.transect = {"autour de La Réunion, Maurice et Tromelin";
      'autour de La Reunion et de Mauritius';
      'autour de Tromelin';
      'autour du sud ouest de La Reunion'};
  case 'brm'
    subset.name = {'img_entiere'; 'Mayotte_Glorieuses'; 'Mayotte_DarEsSalaam'; 'JuanDeNova'};
    subset.limits = [full_img; % #1 img_entiere
      -10.1808 49.1807 -17.4298 42.0793; % #2 Mayotte/Glorieuses
      -13.6117 37.7386 -5.1433 48.0842; % #3 Mayotte/Dar es Salaam
      -20.1347 37.7386 -13.6117 43.8084]; % #4 Juan de Nova
    email_text.sbj = {'Cartes Bougainville - Mayotte';
      'zoom Mayotte - Glorieuses';
      'zoom Mayotte - Dar es Salaam';
      'zoom Juan de Nova'};
    email_text.transect = {'autour de Mayotte, Glorieuses, DarEsSalaam, et Juan-de-Nova';
      'autour de Mayotte et des Glorieuses';
      'le long du transect entre Mayotte et Dar es Salaam';
      'autour de Juan de Nova'};
  case 'bre'
    subset.name = {'img_entiere'; 'JuanDeNova'; 'Europa'; 'SudMada'};
    subset.limits = [full_img; % #1 img_entiere
      -18.1347 41.6812 -16.0757 43.8084; % #2 Juan de Nova
      -23.3354 38.9541 -18.031 42.3357; % #3 Europa
      full_img(1) 39.6086 -22.0048 full_img(4)]; % #4 Sud Mada
    email_text.sbj = {'Cartes Bougainville - Juan de Nova, Europa, et Sud Madagascar';
      'zoom Juan de Nova';
      'zoom Europa';
      'Sud Madagascar'};
    email_text.transect = {'autour de Juan de Nova, Europa, et au Sud de Madagascar';
      "autour de Juan de Nova";
      "autour de Europa";
      "le long du transect retour au sud de Madagascar"};
  case 'bcv'
    subset.name = {'img_entiere'; 'NC_Vanuatu'; 'Chesterfield'};
    subset.limits = [full_img; % #1 img_entiere
      -14.3728 170.3449 -23.5369 162.9219; % #2 NC - Vanuatu
      -17 156 -25 163.5]; % #3 Chesterfield
    email_text.sbj = {'Cartes Bougainville - NC, Vanuatu, et Chesterfield';
      'zoom NC - Vanuatu';
      'zoom Chesterfield';};
    email_text.transect = {"autour de la Nouvelle-Calédonie, des Vanuatu, et Chesterfield";
      "autour de la Nouvelle Calédonie et des Vanuatu";
      "autour de l'archipel de Chesterfield"};
  case 'bcs'
    subset.name = {'img_entiere', 'Vanuatu_Salomon'};
    subset.limits = [full_img; % #1 img_entiere
      -16.4004 159.3371 -7.3472 168.3678]; % #2 Vanuatu - Salomon
    email_text.sbj = {'Cartes Bougainville - Vanuatu et Salomons';
      'zoom Vanuatu - Salomon';};
    email_text.transect = {"autour des Vanuatu et des Iles Salomons"; 
      'autour des Vanuatu et des Iles Salomons';};
  case 'bcp'
    subset.name = {'img_entiere'; 'Salomon_Papouasie'; 'PNG_NC'};
    subset.limits = [full_img; % #1 img_entiere
      -12.1104 144.1371 -5.7772 160.7378; % #2 Salomon - PNG
      -21.8327 150.6626 -10.1355 164.7955]; % #3 PNG - NC
    email_text.sbj = {'Cartes Bougainville - Vanuatu et Salomons';
      'zoom Salomon - PNG';
      'zoom PNG - NC'};
    email_text.transect = {"autour des Vanuatu et des Iles Salomons";
      'autour des Salomon et de la Papouasie Nouvelle Guinée';
      'le long du transect entre la Papouasie Nouvelle Guinée et la Nouvelle Calédonie'};
  case 'bcc'
    subset.name = {'img_entiere'; 'NC_ouest'; 'NC_nord'};
    subset.limits = [full_img; % #1 img_entiere
      -25.3925 160.6854 -19.4450 167.7799; % #2 Ouest
      -23.5 162 -16.5 169]; % #3 Nord
    email_text.sbj = {'Cartes Bougainville - NC';
      'zoom - NC Ouest';
      'zoom NC Nord'};
    email_text.transect = {"autour de la Nouvelle Calédonie"; 
      "à l'Ouest de la Nouvelle Calédonie"; 
      "au Nord de la Nouvelle Calédonie"};
  case 'bcf'
    subset.name = {'img_entiere'; 'Fiji_sudouest'; 'Fiji_nordest'};
    subset.limits = [full_img; % #1 img_entiere
      -25 170 -16 180; % #2 Sud-Ouest
      -19 177 -5 172]; % #3 Nord-Est
    email_text.sbj = {'Cartes Bougainville - Fidji-Tonga';
      'zoom - Fidji Sud-Ouest';
      'zoom Fiji Nord-Est'};
    email_text.transect = {"autour de Fidji-Tonga"; 
      "au Sud-Ouest de Fdji";
      "au Nord-Est de Fidji"};
  otherwise
    subset.name = {'img_entiere'};
    subset.limits = full_img;
    if strcmp(isl_id(2), 'r')
      email_text.sbj = 'Cartes Bougainville - La Reunion';
      email_text.transect = "autour de La Réunion, Mayotte et des Iles Eparses de l'Ocean Indien";
    elseif strcmp(isl_id(2), 'c')
      email_text.sbj = 'Cartes Bougainville - NC';
      email_text.transect = "autour de la Nouvelle Caledonie, des Iles Salomon, de la Papouasie Nouvelle Guinée, ou des Fidji";
    end
end

%% check number of subset to plot
if nargin < 2
  subset_num = 1:size(subset, 1);
  maptype = 'fast';
  variable_toplot = 'chlor_a_poly';
elseif nargin < 3
  if max(subset_num) > size(subset, 1)
    subset_num = 1:size(subset, 1);
  end
  maptype = 'fast';
  variable_toplot = 'chlor_a_poly';
elseif nargin < 4
  if max(subset_num) > size(subset, 1)
    subset_num = 1:size(subset, 1);
  end
  variable_toplot = 'chlor_a_poly';
else
  if max(subset_num) > size(subset, 1)
    subset_num(subset_num > size(subset, 1)) = [];
  end
  if isempty(subset_num)
    error('subset # asked > number of subset defined')
  end
end

%% make sure the lat lon are in the right direction
if abs(mean(diff(lat(1,:)))) > abs(mean(diff(lat(:,1))))
  fprintf('Latitude gradient detected on dim 2 and lon on dimension 1, transposing %s data\n', img_name)
  lat = lat';
  lon = lon';
  for tr = 1:size(data_out,1)
    data_out.data_avg{tr} = data_out.data_avg{tr}';
    data_out.data_std{tr} = data_out.data_std{tr}';
    data_out.data_n{tr} = data_out.data_n{tr}';
  end
end

%% load bathymetry data from GEBCO (full resolution)
if isfile(fullfile(path_containg_img, [isl_id '_bathymetry.mat']))
  fprintf('%s bathymetry data loaded\n', isl_id)
  load(fullfile(path_containg_img,[isl_id '_bathymetry.mat']), 'bathymetry')
else
  fprintf('Extracting %s bathymetry data ... \n', isl_id)
  bathymetry = extractBathymetry(lat, lon, path_bathymetry);
  save(fullfile(path_containg_img,[isl_id '_bathymetry.mat']), 'lat', 'lon', 'bathymetry')
end

mercator_files = {dir(fullfile(path_containg_img,['cmems_mod_glo_*_0.083deg_P1D-m_*m_' isl_id '.nc'])).name}';
% load mercator currents u0 and v0
if any(contains(mercator_files, 'phy-cur_anfc'))
  [~, lat_merc, lon_merc, ~, u_single] = readMERCATOR(fullfile(path_containg_img,...
    mercator_files{contains(mercator_files, 'phy-cur_anfc')}), 'uo', 1, [min(selOC.dt) max(selOC.dt)]);
  % load mercator currents 
  [~, ~, ~, ~, v_single] = readMERCATOR(fullfile(path_containg_img,...
    mercator_files{contains(mercator_files, 'phy-cur_anfc')}), 'vo', 1, [min(selOC.dt) max(selOC.dt)]);
  current_bool = true;
elseif any(contains(mercator_files, 'phy_my_'))
  [~, lat_merc, lon_merc, ~, u_single] = readMERCATOR(fullfile(path_containg_img,...
    mercator_files{contains(mercator_files, 'phy_my_')}), 'uo', 1, [min(selOC.dt) max(selOC.dt)]);
  % load mercator currents 
  [~, ~, ~, ~, v_single] = readMERCATOR(fullfile(path_containg_img,...
    mercator_files{contains(mercator_files, 'phy_my_')}), 'vo', 1, [min(selOC.dt) max(selOC.dt)]);
  current_bool = true;
else
  current_bool = false;
end
% load mercator mixed-layer depth
% if any(contains(mercator_files, 'phy_anfc')) && mld_bool
%   [~, lat_merc, lon_merc360, ~, mld] = readMERCATOR(fullfile(path_containg_img,...
%     mercator_files{contains(mercator_files, 'phy_anfc')}), 'mlotst', 1, [min(selOC.dt) max(selOC.dt)]);
%   mld_bool = true;
% elseif any(contains(mercator_files, 'phy_my_'))
%   [~, lat_merc, lon_merc360, ~, mld] = readMERCATOR(fullfile(path_containg_img,...
%     mercator_files{contains(mercator_files, 'phy_my_')}), 'mlotst', 1, [min(selOC.dt) max(selOC.dt)]);
%   mld_bool = true;
% else
%   mld_bool = false;
% end
mld_bool = false;
lon_merc360 = lon_merc;
lon_merc(lon_merc > 180) = lon_merc(lon_merc > 180)-360;
lon_merc(lon_merc < -180) = lon_merc(lon_merc < -180)+360;

for s = subset_num
  % apply subset: find idx of subset row and col
  idrow = [find(abs(lat(:,1) - subset.limits(s,1)) == min(abs(lat(:,1) - subset.limits(s,1)))); ...
    find(abs(lat(:,1) - subset.limits(s,3)) == min(abs(lat(:,1) - subset.limits(s,3))))];
  idcol = [find(abs(lon(1,:) - subset.limits(s,2)) == min(abs(lon(1,:) - subset.limits(s,2)))); ...
    find(abs(lon(1,:) - subset.limits(s,4)) == min(abs(lon(1,:) - subset.limits(s,4))))];
  if current_bool
    idrow_curr = [find(abs(lat_merc(:,1) - subset.limits(s,1)) == min(abs(lat_merc(:,1) - subset.limits(s,1)))); ...
      find(abs(lat_merc(:,1) - subset.limits(s,3)) == min(abs(lat_merc(:,1) - subset.limits(s,3))))];
    idcol_curr = [find(abs(lon_merc(1,:) - subset.limits(s,2)) == min(abs(lon_merc(1,:) - subset.limits(s,2)))); ...
      find(abs(lon_merc(1,:) - subset.limits(s,4)) == min(abs(lon_merc(1,:) - subset.limits(s,4))))];
  end
  attachements = cell(size(variable_toplot, 2), 2);
  skip_image = false(size(attachements, 1), 1);
  for v = 1:size(variable_toplot, 2)
    id_vartoplot = strcmp(data_out.variable_name, variable_toplot{v});
    prod = data_out.data_avg{id_vartoplot}(min(idrow):max(idrow),min(idcol):max(idcol));
    bathymetrys = bathymetry(min(idrow):max(idrow),min(idcol):max(idcol));
    lats = lat(min(idrow):max(idrow),min(idcol):max(idcol));
    lons = lon(min(idrow):max(idrow),min(idcol):max(idcol));
    if current_bool
      u_single_s = u_single(min(idrow_curr):max(idrow_curr),min(idcol_curr):max(idcol_curr));
      v_single_s = v_single(min(idrow_curr):max(idrow_curr),min(idcol_curr):max(idcol_curr));
      lat_merc_s = lat_merc(min(idrow_curr):max(idrow_curr),min(idcol_curr):max(idcol_curr));
      lon_merc_s = lon_merc360(min(idrow_curr):max(idrow_curr),min(idcol_curr):max(idcol_curr));
      % prepare current data to plot: subsample current data to reduce the number of arrow for m_quiver
      ylight = 1:ceil(size(lat_merc_s,1)/50):size(lat_merc_s, 2);
      xlight = 1:ceil(size(lon_merc_s,2)/50):size(lon_merc_s, 1);
    end
    % set email subject and text
    text_tosend = sprintf("Salut %s,\nVoici les dernières cartes satellites %s " + ...
      "avec des suggestion(s) de station(s).\nBonne mer,\nGuillaume", voab, email_text.transect{s});
  
    % remove shallow pixels:
    prod(bathymetrys >= -30) = NaN;
    
    % wrap station lon
    station.lon360 = station.lon;
    lons360 = lons;
    if all(any(lons(:) < 0) & any(lons(:) >= 0) & ...
        sum(lons(:) < -90 | lons(:) > 90) > sum(lons(:) > -90 & lons(:) < 90))
      station.lon360(station.lon360 < 0) = station.lon(station.lon < 0) + 360;
      lons360(lons360 < 0) = lons360(lons360 < 0) + 360;
    end
    stations = station(station.lat > min(lats(:)) & station.lat < max(lats(:)) & ...
      station.lon360 > min(lons360(:)) & station.lon360 < max(lons360(:)), :);
    
    
    if current_bool || mld_bool

      % ylight = 1:ceil(1/(size(lat_merc_s,1)/size(lats, 1))):size(lat_merc_s, 2);
      % xlight = 1:ceil(1/(size(lat_merc_s,2)/size(lons, 1))):size(lat_merc_s, 1);

      % ylight = 1:ceil(1/(size(lat_merc_s,1)/size(lats, 1)/4)):size(lat_merc_s, 2);
      % xlight = 1:ceil(1/(size(lat_merc_s,2)/size(lons, 1)/4)):size(lat_merc_s, 1);

      ylight = 1:ceil(size(lat_merc_s,1)/40):size(lat_merc_s, 2);
      xlight = 1:ceil(size(lon_merc_s,2)/40):size(lon_merc_s, 1);

      % ylight = 1:ceil(1/size(lats, 1)*100):size(lat_merc_s, 2);
      % xlight = 1:ceil(size(lat_merc_s,2)/size(lons, 1)*5):size(lat_merc_s, 1);

      % ylight = 1:ceil(1*max(size(lat_merc))/size(lats, 1)/4):size(lat_merc, 2);
      % xlight = 1:ceil(1*max(size(lat_merc))/size(lats, 1)/4):size(lat_merc, 1);
    end
    
    % plot
    screen_sz = get(0,'screensize');
    if screen_sz(3) == 2048
      plot_w = 0.3;
      plot_h = 0.5;
    elseif screen_sz(3) == 2560
      plot_w = 0.25;
      plot_h = 0.4;
    elseif screen_sz(3) == 3440
      plot_w = 0.2;
      plot_h = 0.4;
    else
      plot_w = 0.5;
      plot_h = 0.5;
    end
    figh = figure('units','normalized','outerposition',[0 0.025 plot_w plot_h]);
    [~, cbh] = QuickMap(prod, lats, lons, 'viridis', prctile(prod(:), [10 90]), 'log', figh, [], 'f', 'gshhs', maptype); % morgenstemning viridis
    [cbh.Label.String, ~, title_short] = getLabel(variable_toplot{v});
    hold on % hold on to add surface current on top of the map
    
    p = [];
    % load eez data
    eez = shaperead('eez_boundaries_v11.shp'); % [xmin,ymin;xmax,ymax]
    % plot eez borders
    for g = 1:size(eez, 1)
      loneez360 = eez(g).X;
      if all(any(lons(:) < 0) & any(lons(:) >= 0) & ...
          sum(lons(:) < -90 | lons(:) > 90) > sum(lons(:) > -90 & lons(:) < 90))
        loneez360(loneez360 < 0) = loneez360(loneez360 < 0) + 360;
      end
      if strcmp(maptype, 'fast')
        p(1) = plotm(eez(g).Y, loneez360, 'Color', 'r', 'Linewidth', 1);
      elseif strcmp(maptype, 'fancy')
        p(1) = m_plot(loneez360, eez(g).Y, 'Color', 'r', 'Linewidth', 1);
      end
    end
    eez_12m = shaperead('eez_12nm_v3.shp'); % [xmin,ymin;xmax,ymax]
    % plot eez 12m borders
    for g = 1:size(eez_12m, 1)
      lon12meez360 = eez_12m(g).X;
      if all(any(lons(:) < 0) & any(lons(:) >= 0) & ...
          sum(lons(:) < -90 | lons(:) > 90) > sum(lons(:) > -90 & lons(:) < 90))
        lon12meez360(lon12meez360 < 0) = lon12meez360(lon12meez360 < 0) + 360;
      end
      if strcmp(maptype, 'fast')
        p(2) = plotm(eez_12m(g).Y, lon12meez360, 'Color', 'w', 'Linewidth', 1);
      elseif strcmp(maptype, 'fancy')
        p(2) = m_plot(lon12meez360, eez_12m(g).Y, 'Color', 'w', 'Linewidth', 1);
      end
    end
    
    if strcmp(maptype, 'fast')
      geoshow(stations.lat, stations.lon360,'DisplayType','Point','Marker','+','Color','red','Markersize',15,'LineWidth',2);
      textm(stations.lat, stations.lon360, stations.name, 'Color', 'r', 'FontWeight','bold', 'FontSize', 12,'LineWidth',2);
      if current_bool
        quivermc(lat_merc_s(xlight, ylight), lon_merc_s(xlight, ylight), u_single_s(xlight, ylight), ...
          v_single_s(xlight, ylight), 'color', 'w', 'linewidth', 1.8)
        quivermc(lat_merc_s(xlight, ylight), lon_merc_s(xlight, ylight), u_single_s(xlight, ylight), ...
          v_single_s(xlight, ylight), 'color', 'k', 'linewidth', 1)
      end
    elseif strcmp(maptype, 'fancy')
      m_scatter(stations.lon360, stations.lat, 200, 'r', 'Marker','+','LineWidth',2)
      m_text(stations.lon360, stations.lat, stations.name, 'Color', 'r', 'FontWeight','bold', 'FontSize', 12,'LineWidth',2);
      % plot surface currents on top of the map (in black and white contour to see them better)
      if current_bool
        m_quiver(lon_merc_s(xlight, ylight), lat_merc_s(xlight, ylight), u_single_s(xlight, ylight), ...
          v_single_s(xlight, ylight), 'w', 'linewidth', 2)
        m_quiver(lon_merc_s(xlight, ylight), lat_merc_s(xlight, ylight), u_single_s(xlight, ylight), ...
          v_single_s(xlight, ylight), 'k', 'linewidth', 1)
      end
    end
    im_name_split = strsplit(img_name, '_');
    title(sprintf('%s %s %s', im_name_split{1}, im_name_split{2}, im_name_split{3}))
    legend(p, 'EEZ border', '12nm band', 'FontSize', 12, 'Location', 'BestOutside', 'color', 'k', 'TextColor', 'w') % NorthEast SouthEast
    name_rawmap = fullfile(path_containg_img, ['map_' title_short '_' subset.name{s} '_' strrep(img_name, '.mat', '.jpg')]);
    saveGraph(strrep(name_rawmap, '.jpg', ''), 'jpg', 200);
    answer = questdlg('Add datatips on map to highlight sampling locations, press "s" to save and continue','Sampling suggestion','Ok', 'Skip', 'Ok');
    % ask user input to press s to save sampling locations
    if strcmp(answer, 'Ok')
      % populate attachement cellstring and save raw map
      attachements{v,1} = name_rawmap;
      k = '';
      while ~strcmp(k, 's')
        w = waitforbuttonpress();
        if w
          k = get(gcf, 'CurrentCharacter');
          if ~strcmp(k, 's')
            fprintf('?\n');
          end
          if strcmp(k, 's')
            fprintf('Saving map and sending emails ...\n');
          end
        end
      end
      % populate attachement cellstring and save map with sampling station suggestions
      attachements{v,2} = fullfile(path_containg_img, ['map_stations_' title_short '_' subset.name{s} '_' strrep(img_name, '.mat', '.jpg')]);
      saveGraph(strrep(attachements{v,2}, '.jpg', ''), 'jpg', 200);
    elseif strcmp(answer, 'Skip')
      skip_image(v) = true;
      warning('Image skipped')
    end
    close gcf
  end
  % send email to all receipiants
  attachements(skip_image, :) = [];
  attachements = attachements(:);
  if any(~cellfun('isempty', attachements))
    for to = 1:size(emailto,2)
      send_email(emailfrom, pwd, emailto{to}, ['Cartes satellites: ' email_text.sbj{s} ' (email automatique)'], text_tosend, attachements)
      fprintf('Email sent to %s. Subject: %s\n', emailto{to}, ['Cartes satellites: ' email_text.sbj{s} ' (email automatique)'])
    end
  else
    warning('Email not sent. Subject: %s', ['Cartes satellites: ' email_text.sbj{s} ' (email automatique)'])
  end
end
