% Download_Plot_AMISR_from_madrigal.m

% Author: P. Reyes
% Date: July 23, 2025

% What this code does:
% 1) Uses the madrigal API to fetch and list the available instruments that
% match the name 'IS Radar' and shows the instrument code, name, latitude,
% and longitude.
% 2) Selects the PFISR instrument and finds the experiments that
% happend in the specified time range, and list their names, description,
% url, and experiment id
% 3) List the files associated with the selected experiment
% 4) Selects and download a a specific file by searching a kindat
% description that matches e.g. Long Pulse and 5 min
% 5) Reads the downloade file and display information about the experiment 
% 7) Display information about the arrays
% 8) Plots the electron density of all the beams.

% This file uses the following files inside a folder named madmatlab
% The required files can be downloaded from the cedar madrigal repository:
% https://cedar.openmadrigal.org/docs/name/rt_matlab.html

% required files:

% getMadrigalCgiUrl.m
% getInstrumentsWeb.m
% getExperimentsWeb.m
% getExperimentFilesWeb.m
% madDownloadFile.m

addpath(path('madmatlab',path)); % to find the required files


madurl = 'https://data.amisr.com/madrigal/';

cgiurl = getMadrigalCgiUrl(madurl);

fprintf('List all IS Radars, and their latitudes and longitudes:\n');
instArray = getInstrumentsWeb(cgiurl);
for i = 1:length(instArray)
    if contains(instArray(i).name,'IS Radar','IgnoreCase',true)
        fprintf('Instr. code: %i, instr.name: %s, at lat %f and long %f\n', ...
            instArray(i).code, ...
            instArray(i).name, ...
            instArray(i).latitude, ...
            instArray(i).longitude);
    end
end

PFISR_code = 61;

% now list all experiments from local Madrigal site with PFISR (code 61) in
% 30 Aug, 2014 - should be data if default files installed...
fprintf('Find experiments in a time range:\n');

startdate = datenum('7/22/2025');
enddate = datenum('7/23/2025');
expArray = getExperimentsWeb(cgiurl, PFISR_code, startdate, enddate, 1);
for i = 1:length(expArray)
    fprintf('Experiment name: %s, at url %s and id %i\n', ...
            expArray(i).name, ...
            expArray(i).url, ...
            expArray(i).id);
end

% now list all files in the first experiment
fprintf('List all files in the first experiment (IPY27):\n');

% selecting 
expFileArray = getExperimentFilesWeb(cgiurl, expArray(1).id);
for i = 1:length(expFileArray)
    fprintf('%i:File name: %s, kindat %i, descr: %s\n', i,...
            expFileArray(i).name, ...
            expFileArray(i).kindat,...
            expFileArray(i).kindatdesc);
    if (contains(expFileArray(i).kindatdesc,'Fitted - Long Pulse', ...
            IgnoreCase=true) & ... 
            contains(expFileArray(i).kindatdesc,'5 min', ...
            IgnoreCase=true))
        lp_file = expFileArray(i).name;
    end
end

saved_hdf5file = 'pfa_IPY27_22Jul2025_1.h5';
% download the Long Pulse data
fprintf('\nDownload file : %s\n',lp_file);

result = madDownloadFile(cgiurl, lp_file, ...
    saved_hdf5file, ...
    'ISR School', 'isr.summer.school@gmail.com', ...
         'ISR Summer School 2025', 'hdf5');


% read HDF5 file
fprintf('\nDisplaying HDF5 file information.\n')
h5disp(saved_hdf5file,'/','min') % 'min', displays only the group and dataset names

h5disp(saved_hdf5file,'/Metadata','min');

PFISR_expinfo = h5read(saved_hdf5file,'/Metadata/Experiment Parameters')

[nrows, ncols] = size(PFISR_expinfo.name);

fprintf('\nDisplay information about the Experiment:\n\n');
for i=1:nrows
    fprintf('%s : %s\n',char(PFISR_expinfo.name(i,:)), ...
        char(PFISR_expinfo.value(i,:))); 
end


fprintf('\nDisplay /Data/Array info with h5disp:\n\n');

h5disp(saved_hdf5file,'/Data/Array Layout','min');


% read DataArray
DataArray = h5info(saved_hdf5file, '/Data/Array Layout');

fprintf('\nDisplay /Data/Array info:\n\n');

keySet = {};
valueSet = {};
allbeams = [];
for i = 1:length(DataArray.Groups)
    fprintf('%i : %s\n',i,DataArray.Groups(i).Name);
    
    for j = 1:length(DataArray.Groups(i).Datasets)
        toread = sprintf("%s/%s",DataArray.Groups(i).Name,DataArray.Groups(i).Datasets(j).Name);
        array = h5read(saved_hdf5file, toread);
        beamname = DataArray.Groups(i).Name(end-5:end-1);
        keySet = [keySet ; [beamname,' ' , DataArray.Groups(i).Datasets(j).Name]];
        valueSet = [valueSet ; array];
        allbeams = [allbeams;beamname];
    end
    for j = 1:length(DataArray.Groups(i).Groups)
        fprintf('    %i: %s\n',j,DataArray.Groups(i).Groups(j).Name);
        for k = 1:length(DataArray.Groups(i).Groups(j).Datasets)
            toread = sprintf("%s/%s",DataArray.Groups(i).Groups(j).Name,DataArray.Groups(i).Groups(j).Datasets(k).Name);
            array = h5read(saved_hdf5file, toread);
            keySet = [keySet ; [DataArray.Groups(i).Name(end-5:end) , DataArray.Groups(i).Groups(j).Datasets(k).Name]];
            valueSet = [valueSet ; array];
            fprintf('      %s\n',DataArray.Groups(i).Groups(j).Datasets(k).Name);
        end
    end
end
PFISR_data = containers.Map(keySet,valueSet);

allbeams = unique(allbeams,'rows');

[nbeams,beamdigits] = size(allbeams);

fprintf('\nPlot Electron density:\n\n');



for row = 1:nbeams
    beam = allbeams(row,:);
    Xtimestamps = double(PFISR_data([beam,' timestamps']));
    Mtimestamps = zeros('like',Xtimestamps); 
    for i=1:length(Xtimestamps),
        Mtimestamps(i) = datenum([1970 1 1 0 0 Xtimestamps(i)]);
    end

    subplot(nbeams,1,row)
    h = pcolor(Mtimestamps,PFISR_data([beam,' range'])/1e3,...
        transpose(PFISR_data([beam,' ne'])));
    clim([1e10,1e12]);
    set(gca,'ColorScale','log')
    datetick('x','keeplimits');
    set(h, 'EdgeColor', 'none');
    colorbar;
    azm = PFISR_data([beam,' azm']);
    elm = PFISR_data([beam,' elm']);
    title(sprintf('PFISR Ne, %s (%.2f az, %.2f el)\n%s - %s',...
        beam,azm(1), elm(1),...
        datestr(Mtimestamps(1)),datestr(Mtimestamps(end)) ))
    xlabel("Time (UT)")
    ylabel("Range (km)")
end


saveas(gcf,'PFISR_IPI27_22Jul2025_1.png');