%% Elaborate raw input files from the Broadband Radiation Probe @MOSAiC
%% and export to Pangaea terrestrial and solar files.
%
% Michael Lonardi
% Contact: michael.lonardi@uni-leipzig.de
%
%% SELECT DAY (uncomment needed flight)
%flight_day = 20200629.1; % flight 1
%flight_day = 20200629.2; % flight 2
%flight_day = 20200710.4; % flight 4
%flight_day = 20200713.6; % flight 6
%flight_day = 20200714.8; % flight 8
%flight_day = 20200715.9; % flight 9
%flight_day = 20200715.10; % flight 10
%flight_day = 20200715.11; % flight 11
%flight_day = 20200719.14; % flight 14
%flight_day = 20200722.20; % flight 20 
%flight_day = 20200723.22; % flight 22
%flight_day = 20200723.24; % flight 24
%flight_day = 20200724.27; % flight 27
%flight_day = 20200725.30; % flight 30
%flight_day = 20200726.31; % flight 31 and 32
%flight_day = 20200726.32; % flight 32
%flight_day = 20200727.331; % flight 33 first part
%flight_day = 20200727.332; % flight 33 second part

switch flight_day
    case 20200629.1 % June 29 morning === flight 1
        path.date = 20200629;
        path.dateshort = 200629;
        path.time = 075332;
        path.timeshort = 0753;
        flight_number = 1;
        startingOffset = 10;
        final = 83489; % end - 10
        % load('200629-f1'); % final dataset 
        start_ascent = 35728; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 63669;
    case 20200629.2 % June 29 afternoon === flight 2 
        path.date = 20200629;
        path.dateshort = 200629;
        path.time = 121515;
        path.timeshort = 1215;
        flight_number = 2;
        startingOffset = 10;
        final = 89096; % end -10
        % load('200629-f2'); % final dataset
        start_ascent = 24547; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 71063;
    case 20200710.4 % July 10 === flight 4
        path.date = 20200710;
        path.dateshort = 200710;
        path.time = 120200;
        path.timeshort = 1202;
        flight_number = 4;
        startingOffset = 10;
        final = 36855; % end -10
        % load('200710-f4'); % final dataset
        start_ascent = 8962; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 29518;
    case 20200713.6 % July 13 === flight 6
        path.date = 20200713;
        path.dateshort = 200713;
        path.time = 133046;
        path.timeshort = 1330;
        startingOffset = 20; 
        final = 84194;
        flight_number = 6;
        % load('200713-f6'); % final dataset
        start_ascent = 6701; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 78330;
    case 20200714.8 % July 14 === flight 8 === SINGLE LAYER CLOUD starting from the surface
        path.date = 20200714;
        path.dateshort = 200714;
        path.time = 132046;
        path.timeshort = 1320;
        startingOffset = 700; % also in-flight problems with height calculation, run manually
        final = 74223; % end - 10
        flight_number = 8;
        % load('200714-f8'); % final dataset
        start_ascent = 6204; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 60680;
    case 20200715.9 % July 15 === flight 9
        path.date = 20200715;
        path.dateshort = 200715;
        path.time = 073516;
        path.timeshort = 0735;
        startingOffset = 100; 
        final = 71100; % end - 10;
        flight_number = 9;
        % load('200715-f9'); % final dataset
        start_ascent = 11568; % starts from 10 m, gaps in the first 60 m
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 55517; % ends at 15 m because the end is a long slope
    case 20200715.10 % July 15 === flight 10 alone
        path.date = 20200715;
        path.dateshort = 200715;
        path.time = 112544;
        path.timeshort = 1125;
        startingOffset = 1000; 
        final = 63388; % 
        flight_number = 10;
        % load('200715-f10'); % final dataset
        start_ascent = 19120; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 63388; 
    case 20200715.11 % July 15 === flight 11 alone, careful with indexes
        startingOffset = 63388 +1; 
        final = 76742; % 
        flight_number = 11;
        % load('200715-f11'); % final dataset
        start_ascent = 63388 +1; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 76742; 
    case 20200719.14 % July 19 === flight 14
        path.date = 20200719;
        path.dateshort = 200719;
        path.time = 130036;
        path.timeshort = 1300;
        startingOffset = 22;
        final = 68917; % end -10
        flight_number = 14;
        % load('200719-f14'); % final dataset
        start_ascent = 6540; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 63062;   
     case 20200722.20 % July 22 === flight 20
        path.date = 20200722;
        path.dateshort = 200722;
        path.time = 033922;
        path.timeshort = 0339;
        startingOffset = 10;
        final = 73392;
        flight_number = 20;
        % load('200722-f20'); % final dataset
        start_ascent = 3460; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 62158; 
    case 20200723.22 % July 23 morning === flight 22
        path.date = 20200723;
        path.dateshort = 200723;
        path.time = 074509;
        path.timeshort = 0745;
        startingOffset = 10;
        final = 53273; % end -10
        flight_number = 22;
        % load('200723-f22'); % final dataset
        start_ascent = 4454; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 49061; 
    case 20200723.24 % July 23 afternoon === flight 24 
        path.date = 20200723;
        path.dateshort = 200723;
        path.time = 110032;
        path.timeshort = 1100;
        startingOffset = 10; 
        final = 74693; % end-10
        flight_number = 24;
        % load('200723-f24'); % final dataset
        start_ascent = 8536; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 67235; 
    case 20200724.27 % July 24 morning === flight 27
        path.date = 20200724;
        path.dateshort = 200724;
        path.time = 085403;
        path.timeshort = 0854; 
        startingOffset = 10;  
        final = 96113; % end-10
        flight_number = 27;
        % load('200724-f27'); % final dataset
        start_ascent = 9230; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 94589; 
    case 20200725.30 % July 25 === flight 30
        path.date = 20200725;
        path.dateshort = 200725;
        path.time = 134656;
        path.timeshort = 1346; 
        startingOffset = 10; % guess
        final = 63076; % end-10
        flight_number = 30;
        % load('200725-f30') % final dataset
        start_ascent = 13650; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 56779; 
    case 20200726.31 % July 26 === flight 31 only
        path.date = 20200726;
        path.dateshort = 200726;
        path.time = 112854;
        path.timeshort = 1128;
        startingOffset = 10; 
        final = 91419;
        flight_number = 31;
        % load('200726-f32') % final dataset
        start_ascent = 25528; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 91419; 
    case 20200726.32 % July 26 === flight 32 only
        startingOffset = 91419 +1; 
        final = 161514;
        flight_number = 32;
        % load('200726-f32') % final dataset
        start_ascent = 100743; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 156291;
    case 20200727.331 % July 27 === flight 33 first part
        path.date = 20200727;
        path.dateshort = 200727;
        path.time = 072823;
        path.timeshort = 0728;
        startingOffset = 10; 
        final = 67554;
        flight_number = 33;
        % load('200727-f33a') % final dataset
        start_ascent = 10339; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 67554;
    case 20200727.332 % July 27 === flight 33 second part
        startingOffset = 67554 +1; 
        final = 81676;
        flight_number = 34;
        % load('200727-f33b') % final dataset
        start_ascent = 67554 +1; 
        stop_ascent = hmax_index;
        start_descent = hmax_index +1;
        stop_descent = 76288;
end

% startingOffset and final are used to skip the incomplete lines at the
% beginning and end of the files.
% ascent/descent are used to crop the meaninful part of the flights, since
% the instrument often had an idle time on the ground before and after the
% flight

%% Create paths 
if path.timeshort > 959 % (later than 9:59) from 10 AM everyting is fine
    path.bme = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\BME_1_',num2str(path.date),'_',num2str(path.time),'.txt']);
    path.ads = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\ads1115_1_',num2str(path.date),'_',num2str(path.time),'.txt']);
    path.adc = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\',num2str(path.dateshort),'_',num2str(path.timeshort),'_Data_ADC.txt']);
    path.gps = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\GPS_1_',num2str(path.date),'_',num2str(path.time),'.txt']);
    path.dof = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\dof10_1_',num2str(path.date),'_',num2str(path.time),'.txt']);
else % if time is less than 10 AM, the 0 in the time (e.g. 08:59 --> 8:59) went lost and we have to add it manually
    path.bme = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\BME_1_',num2str(path.date),'_0',num2str(path.time),'.txt']);
    path.ads = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\ads1115_1_',num2str(path.date),'_0',num2str(path.time),'.txt']);
    path.adc = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\',num2str(path.dateshort),'_0',num2str(path.timeshort),'_Data_ADC.txt']);
    path.gps = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\GPS_1_',num2str(path.date),'_0',num2str(path.time),'.txt']);
    path.dof = (['C:\Users\Michael\Desktop\MOSAiC_Data\MOSAiC_Radiation\BP1\',num2str(path.date),'\dof10_1_',num2str(path.date),'_0',num2str(path.time),'.txt']);
end

%% Read the raw data from the text files
DataBME = importBME(path.bme,1,inf); % T/P sensors
DataADC = importADC(path.adc,1,inf);
DataADS = importADS(path.ads,1,inf);
DataDOF = importDOF(path.dof,[1, Inf]); % attitude data
DataGPS = importGPS(path.gps,[1, Inf]); % GPS data

%% CONSTANTS
% Atmospheric constants
c_p  = 1005.0; % specific heat in J kg-1 K-1
dens = 1.2922; %kg m-3 for 0C Density of air 
bolzmann = 5.67e-8;

% ads1115 constants
GAIN = 1;
data_rate = 128;
r0 = 40000; % ohm

% conversion factor counts to voltage
vocofac = (2*4.096/GAIN)/(32767--32768);

% sensivity
sens = [10.07,10.25,10.72,10.37]; % pyrano pyrano pyrgeo pyrgeo of bp1 (from calibration papers)
sens = sens*10^(-6); % in V / (W/m²)

%% Times
% Create times vectors starting from the sensors' time and the ADS date
ADCtime = datetime(DataADS.yyyy(1),DataADS.month(1),DataADS.dd(1),DataADC.hh,DataADC.mm,0,DataADC.sssss*1000); 
ADStime = datetime(DataADS.yyyy(1),DataADS.month(1),DataADS.dd(1),DataADS.hh,DataADS.mm,0,DataADS.sssss*1000); 
BMEtime = datetime(DataBME.yyyy(1),DataBME.month(1),DataBME.dd(1),DataBME.hh,DataBME.mm,0,DataBME.sssss*1000); 
DOFtime = datetime(DataDOF.yyyy(1),DataDOF.month(1),DataDOF.dd(1),DataDOF.hh,DataDOF.mm,0,DataDOF.sssss*1000);
GPStime = datetime(DataGPS.yyyy(1),DataGPS.month(1),DataGPS.dd(1),DataGPS.hh,DataGPS.mm,0,DataGPS.sssss*1000);

%% TEMPERATURE CORRECTION from ADS
% Converts the counts from the ADS to voltages
DataADS.u_u = DataADS.digital_u_u .* DataADS.conversion_factor; % voltage of the upward irr
DataADS.u_d = DataADS.digital_u_d .* DataADS.conversion_factor; % voltage of the downward irr
DataADS.u_ref = DataADS.digital_u_ref .* DataADS.conversion_factor; % reference voltage
% Calculate the temperature correction for the pyrgeo body from the ADS file
% Uses the v and v_ref to calculate the thermistor resistance, then from
% that the temperature is calculated.
temp_lw_u = RtoT( VtoR( DataADS.u_u,DataADS.u_ref)); % Temperature of the body for the upward pyrg / d_adc 19
temp_lw_d = RtoT( VtoR( DataADS.u_d,DataADS.u_ref)); % Temperature of the body for the downward pyrg / d_adc 18
% then the temperature is interpolated to the ADC time, as the radiometers
% have a differet time array. The new temperatures are on the same time
% array as the radiometers.
temp_lw_u_interpolated = interp1(ADStime,temp_lw_u,ADCtime); % yy = interp1(x,y,xx)
temp_lw_d_interpolated = interp1(ADStime,temp_lw_d,ADCtime); % yy = interp1(x,y,xx)

%% IRRADIANCES from ADC with T correction from ADS
% Converts from the voltages of ADC to irradiances using the pyranometer
% formula and the sensibility (from calibration) of the respective
% pyranometer. The gain is considered too
sw_d = U_to_SolIrr(DataADC.vswd,sens(1),GAIN); % Downwards solar  irradiance / d_adc 13 
sw_u = U_to_SolIrr(DataADC.vswu,sens(2),GAIN); % Upwards solar  irradiance / d_adc 14 

% for the pyrgeometers applies also the body temperature correction using
% the temperatures from ADS
lw_d = U_to_TerrIrr(DataADC.vlwd,sens(3),GAIN,temp_lw_d_interpolated);
% Downwards terrestrial irradiance / d_adc 16 check

lw_u = U_to_TerrIrr(DataADC.vlwu,sens(4),GAIN,temp_lw_u_interpolated);
% Upwards terrestrial irradiance / d_adc 18 check

% Using the above quantities, calculates the net irradiances
sw_net = sw_d - sw_u; % solar  net irradiance / d_adc 21 
lw_net = lw_d - lw_u; % terrestrial net irradiance / d_adc 22 
tot_net = sw_net + lw_net; % Total net irradiance / d_adc 23 

%% BME
% use the Barometric formula:  h = (T0 / L0) * ((p / P0)**(-(R* * L0) / (g0 * M)) - 1)
% with input in Pa (e.g. 101325 for the surface)
% where:
%       h  = height above sea level
%       T0 = standard temperature at sea level = 288.15
%       L0 = standard temperatur elapse rate = -0.0065
%       p  = measured pressure
%       P0 = static pressure = 1013.25
%       g0 = gravitational acceleration = 9.80665
%       M  = mloecular mass of earth's air = 0.0289644
%       R* = universal gas constant = 8.31432
%       
%       Given the constants, this works out to:
height = (DataBME.temperature(start_ascent)/0.0065) .* (1 - (DataBME.pressure / max(DataBME.pressure(start_ascent:stop_descent))).^0.190263);


% Interpolates the bme data to match with the time array of the radiometers.
temperature_interpolated = interp1(BMEtime,DataBME.temperature,ADCtime); % yy = interp1(x,y,xx)
rh_interpolated = interp1(BMEtime,DataBME.rh,ADCtime); % yy = interp1(x,y,xx)
pressure_interpolated = interp1(BMEtime,DataBME.pressure,ADCtime); % yy = interp1(x,y,xx)
height_interpolated = interp1(BMEtime,height,ADCtime); % yy = interp1(x,y,xx)
% If errors appear here there is one row too much in DataBME and/or in
% BMEtime. Remove it manually and re-run

%% DOF for attitude
% Interpolates the dof data to match with the time array of the radiometers.
roll_interpolated = interp1(DOFtime,DataDOF.roll,ADCtime); % yy = interp1(x,y,xx)
pitch_interpolated = interp1(DOFtime,DataDOF.pitch,ADCtime); % yy = interp1(x,y,xx)
yaw_interpolated = interp1(DOFtime,DataDOF.yaw,ADCtime); % yy = interp1(x,y,xx)

%% Manually check icing and define its array

%% Manually save dataset

%% For the terrestrial radiation, run offline the IDL deconvolution algorithm by Ehrlich et al. (2015) 10.5194/AMT-8-3671-2015 

%% Section: Pangaea export of flights

%% Publishing Terrestrial for PANGAEA (3s running mean)
% Variables to export
Flight.time = ADCtime(start_ascent:stop_descent); % now it is corrected
Flight.lwd = movmean(lw_d_dec(start_ascent-startingOffset:stop_descent-startingOffset),36);
Flight.lwu = movmean(lw_u_dec(start_ascent-startingOffset:stop_descent-startingOffset),36);
Flight.pressure = movmean(pressure_interpolated(start_ascent:stop_descent),36)/100; % in hPa
Flight.temperature = movmean(temperature_interpolated(start_ascent:stop_descent),36)  + 273.15; % in Kelvin
Flight.height = movmean(height_interpolated,36); % Geopotential height2
Flight.rh = movmean(rh_interpolated(start_ascent:stop_descent),36);
Flight.icing = ceil(movmean(icing_flag(start_ascent:stop_descent),36)); % if one over 36 is iced, everything is iced
Flight.time.Format = 'yyyy-MM-dd''T''hh:mm:ss.SSS'; % change of format

% Creates the file
fid = fopen(['PANGAEA/Flight ',num2str(flight_number),' terrestrial radiation.dat'],'w'); % Creates the file
% Header
fprintf(fid, 'DateTime[UTC] F_D[Wm^(-2)] F_U[Wm^(-2)] h[m] p[hPa] T[K] rh[%%] icing[flag]\r\n');
% Body
for i=1:length(Flight.time)
    fprintf(fid, '%s %.6f %.6f %.4f %.2f %.4f %.4f %.0f \r\n', Flight.time(i), Flight.lwd(i), Flight.lwu(i), Flight.height(i), Flight.pressure(i), Flight.temperature(i), Flight.rh(i), Flight.icing(i));
 end
fclose(fid); % Done

%% Publishing Solar for PANGAEA (no running mean)
% Variables to export
Flight.height = height_interpolated2(start_ascent:stop_descent); % Geopotential height
Flight.time = ADCtime(start_ascent:stop_descent); % same as lw now
Flight.swd = sw_d(start_ascent:stop_descent); % raw
Flight.swu = sw_u(start_ascent:stop_descent); % raw
Flight.lat = lat_interpolated(start_ascent:stop_descent); 
Flight.lon = lon_interpolated(start_ascent:stop_descent); 
Flight.roll = roll_interpolated(start_ascent:stop_descent); 
Flight.pitch =pitch_interpolated(start_ascent:stop_descent);
Flight.yaw =yaw_interpolated(start_ascent:stop_descent);
Flight.pressure = pressure_interpolated(start_ascent:stop_descent)/100; % in hPa
Flight.icing = icing_flag(start_ascent:stop_descent); % 
Flight.time.Format = 'yyyy-MM-dd''T''hh:mm:ss.SSS'; % change of format

% Creates the file
fid = fopen(['PANGAEA/Flight ',num2str(flight_number),' solar radiation.dat'],'w'); % Creates the file
% Header
fprintf(fid, 'DateTime[UTC] F_D[Wm^(-2)] F_U[Wm^(-2)] lat[deg] lon[deg] roll[deg] pitch[deg] yaw[deg] h[m] p[hPa] icing[flag]\r\n');
% Body
for i=1:length(Flight.time)
    fprintf(fid, '%s %.6f %.6f %.5f %.5f %.4f %.4f %.4f %.4f %.2f %.0f \r\n', Flight.time(i), Flight.swd(i), Flight.swu(i), Flight.lat(i), Flight.lon(i), Flight.pitch(i), Flight.roll(i), Flight.yaw(i), Flight.height(i), Flight.pressure(i), Flight.icing(i));
 end
fclose(fid); % Done
