function ecuador_flux_main

% Directory where stored
% Use / in path not \
% Also dont forget the trailing '/'
% E.g: H:/xmosaic/tmp/sensor/EcuadorMay10/Matlab/
DIR = 'H:/xmosaic/tmp/sensor/EcuadorMay10/Data/CsvGrid/';

% Read master.csv and compute flux
% numbers for locations in Ecuador
master_fname = 'master.csv';
master_fname = strcat(DIR,master_fname);
master_fmt = '%[^,],%[^,],%[^,],%d,%[^,],%d,%f,%f';

[master.loc master.co2 master.stm1 master.stm1_ind master.stm2 master.stm2_ind master.texture master.bulk] = ...
    textread(master_fname, master_fmt);

for loc=1:length(master.co2),

    % specify file names
    fname_co2 = master.co2(loc);
    fname_stm1 = master.stm1(loc);
    fname_stm2 = master.stm2(loc);

    %fmt of the data
    fmt = '%[^,] %f %f %f %f %f %f %f %f %f %f';

    % read co2 input file
    [co2.dt co2.ch1 co2.ch2 co2.ch3 co2.ch4 co2.ch5 co2.ch6 co2.ch7 co2.ch8 co2.ch9 co2.ch10] = ...
         textread(fname_co2,fmt,'delimiter',',');
    [stm1.dt stm1.ch1 stm1.ch2 stm1.ch3 stm1.ch4 stm1.ch5 stm1.ch6 stm1.ch7 stm1.ch8 stm1.ch9 stm1.ch10] = ...
         textread(fname_stm1,fmt,'delimiter',',');
    [stm2.dt stm2.ch1 stm2.ch2 stm2.ch3 stm2.ch4 stm2.ch5 stm2.ch6 stm2.ch7 stm2.ch8 stm2.ch9 stm2.ch10] = ...
         textread(fname_stm2,fmt,'delimiter',',');

    % convert date string to a number 
    % so that we can compare values
    co2_num =  datetime2datenum(co2.dt, 'yyyy-mm-dd HH:MM:SS');
    stm1_num = datetime2datenum(stm1.dt, 'yyyy-mm-dd HH:MM:SS');
    stm2_num = datetime2datenum(stm2.dt, 'yyyy-mm-dd HH:MM:SS');

    % find the overlap between these three files
    t1 = max([min(co2_num) min(stm1_num) min(stm2_num)]);
    t2 = min([max(co2_num) max(stm1_num) max(stm2_num)]);

    % find the indexes which lie between t1 and t2
    ind_co2 = find(co2_num > t1 & co2_num < t2);
    ind_stm1 = find(stm1_num > t1 & stm1_num < t2);
    ind_stm2 = find(stm2_num > t1 & stm2_num < t2);

    % concatenate Date, CO2, Soil Temp, Soil Moisture
    dt = co2.dt(ind_co2);
    all_modes = [co2.ch1(ind_co2) co2.ch2(ind_co2) co2.ch3(ind_co2) ...
        stm1.ch1(ind_stm1) stm1.ch2(ind_stm1) stm2.ch1(ind_stm2) ...
        stm1.ch3(ind_stm1) stm1.ch4(ind_stm1) stm1.ch3(ind_stm2)];

    % correct the co2 for each depth
    NUM_DEPTH = 3;
    PRESSURE = 1013.25 * 1.0E-5;
    KELVIN_CORRECTION = 273.15;
    GAS_CONSTANT = 8.314;

    co2_corrected = [];
    for i=1:NUM_DEPTH,
        denom = KELVIN_CORRECTION + all_modes(:,i+6);
        co2_corrected(:,i) = all_modes(:,i)*PRESSURE ./ denom;
    end

    [TIMES DEPTHS] = size(co2_corrected);
    GRADIENT = zeros(TIMES,1);

    x = [2, 8, 16];
    CO2_DEPTHS = 1:3;
    for i=1:TIMES
        y = co2_corrected(i,CO2_DEPTHS);    
        p = polyfit(x,y,2);
        k = polyder(p);
        GRADIENT(i) = polyval(k,0);
    end

    %% PLOT THE GRADIENTS
    figure; 
    subplot(2,1,1);
    plot([1:TIMES],GRADIENT);
    tk = floor([1:TIMES/8:TIMES]);
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk));
    title('Flux values');
    subplot(2,1,2); MOISTURE_IND = 5;
    plot([1:TIMES],all_modes(:,MOISTURE_IND));
    tk = floor([1:TIMES/8:TIMES]);
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk));
    title('Moisture');
    ylim([0.3 0.4]);

    DAO = 1.47 * 1.0E-5;
    T_REF = 293.15;
    P_RATIO = 1.028926;

    % Compute Do for temperature values
    TEMP_IND = 8;   % take the depth 8cm temperature
    data_do = DAO*(all_modes(:,TEMP_IND) / T_REF).^1.75 * P_RATIO;
    
    
end


end