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 = 'C:/Ecuador/CsvGrid/';

constants;

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

[master.setup master.co2 master.stm1 master.mind master.stm2 master.sind master.texture master.bulk] = ...
    textread(MASTER_FNAME, MASTER_FMT);

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

    % create filenames
    fname_co2 = strcat(DIR,master.co2(loc),'.csv');
    fname_stm1 = strcat(DIR,master.stm1(loc),'.csv');

    %fmt in which data is stored
    fmt = '%[^,] %f %f %f %f %f %f %f %f %f %f';

    % read co2 input file and then the supporting temp, moisture files
    [co2.dt co2.ch1 co2.ch2 co2.ch3 co2.ch4 co2.ch5 co2.ch6 co2.ch7 co2.ch8 co2.ch9 co2.ch10] = ...
         textread(char(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(char(fname_stm1),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');

    % find the overlap between these three files
    t1 = max([min(co2_num) min(stm1_num)]);
    t2 = min([max(co2_num) max(stm1_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);

    % if the numberof measurements are not the same .. skip
    if (length(ind_co2) ~= length(ind_stm1)) continue; end
    
    % 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) stm1.ch3(ind_stm1) stm1.ch4(ind_stm1)];

    % correct the co2 for each depth
    NUM_DEPTH = 3;
    PRESSURE = 1.01325 * 1.0E5;
    KELVIN_CORRECTION = 273.15;
    GAS_CONSTANT = 8.314;
    TEMP_IND = 3+master.sind(loc);   % take the depth 8cm temperature
    
    co2_corrected = [];
    for i=1:NUM_DEPTH,
        denom = KELVIN_CORRECTION + all_modes(:,TEMP_IND);
        co2_corrected(:,i) = all_modes(:,i)/GAS_CONSTANT*PRESSURE ./ denom;
    end

    % HALF the co2 concentration
    co2_corrected = co2_corrected/2;
    
    [TIMES DEPTHS] = size(co2_corrected);
    gradient = zeros(TIMES,1);
    gradient_lijun = zeros(TIMES,1);
    
    x = [-0.02, -0.08, -0.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);
        % gradient_lijun(i,:)=2*p(i,1)*x+p(i,2);
    end

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

    % Compute Do for temperature values
    data_do = DAO*(( KELVIN_CORRECTION+all_modes(:,TEMP_IND)) / T_REF).^1.75 * P_RATIO;
    data_do_lijun = DAO*((273.15+all_modes(:,TEMP_IND))./(273.15+20)).^1.5;
    
    % Compute the bulk density
    POROSITY_CONSTANT = 2.65;   % the denominator term
    phi = 1-master.bulk(loc)/POROSITY_CONSTANT;
    
    % compute the epsilon for each time step
    MOISTURE_IND = 3+master.mind(loc);
    epsilon = phi - all_modes(:,MOISTURE_IND);
    
    % Millington Quick Model
    dgs_mq = (epsilon.^(10/3).*data_do)/(phi^2);
    % MM model
    pwr = 2.9*master.texture(loc);
    dgs_mm = phi^2.*data_do.*(epsilon/phi).^pwr;
    
    fx_mq = -dgs_mq.*gradient;
    fx_mm = -dgs_mm.*gradient;
    
    %% PLOT THE FLUX
    
    
    % CONC
    figure('position',[100 100 1400 800]);
    subplot(4,1,1); set(gca,'FontSize',font_size,'box','off');
    plot([1:TIMES],co2_corrected);
    tk = floor([1:TIMES/6:TIMES]);
    legend({'2 cm','8 cm','16 cm'});
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk));
    ylabel('Conc [ppm]');
    title(strcat('Concentration -',master.setup(loc)));    
    % MQ
    %{
    subplot(4,1,2); set(gca,'FontSize',font_size,'box','off');
    plot([1:TIMES],fx_mq);
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk)); 
    %ylim([-50 100]);
    title(strcat('Millington Quick -',master.setup(loc)));
    %}
    %MM
    subplot(4,1,2); set(gca,'FontSize',font_size,'box','off');
    plot([1:TIMES],fx_mm);
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk));
    ylim([0 50]);
    title(strcat('Moldrup Model -',master.setup(loc)));
    
    % PLOT THE SOIL MOISTURE
    subplot(4,1,3); set(gca,'FontSize',font_size,'box','off');
    mdata = all_modes(:,MOISTURE_IND);
    ind = find(isnan(mdata) == 1);
    mdata(ind) = 1;
    plot([1:TIMES],mdata);
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk));
    ylabel('Soil Moisture [0..1]');
    sm_min = min([0 ; all_modes(:,MOISTURE_IND)]);
    sm_max = max([0.5 ; all_modes(:,MOISTURE_IND)]);
    if (sm_min < 0.3) sm_min = 0.3; end
    if (sm_max > 0.45) sm_max = 0.45; end
    ylim([sm_min sm_max]);
 
    % PLOT THE SOIL TEMPERATURE
    subplot(4,1,4); set(gca,'FontSize',font_size,'box','off');
    plot([1:TIMES],all_modes(:,TEMP_IND));
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk));
    ylabel('Soil Temperature [celcius]');
    ylim([22 27.5]);
    
    % save figure
    figname = strcat(fig_save_loc, 'flux', master.setup(loc), fig_ext);
    save_figure(gcf, char(figname), fig_fmt, [fig_width fig_height],'points'); 
    
    % Lijuns's code
    %{
    CO2 = all_modes(:,1:3);
    MOIST = all_modes(:,MOISTURE_IND); MOIST = repmat(MOIST,1,3);
    TEMP = all_modes(:,TEMP_IND); TEMP = repmat(TEMP,1,3)
    moldrup = lijun_flux_code(CO2,MOIST,TEMP);
    
    subplot(4,1,4); set(gca,'FontSize',font_size,'box','off');
    plot([1:TIMES],moldrup);
    set(gca,'xtick',tk);
    set(gca,'xticklabel',dt(tk));
    ylabel('Lijun Flux');    
    %}
    

%{
    % output to csv file
    % corrected co2, MQ, MM, moisture and temperature
    out = [co2_corrected fx_mq fx_mm all_modes(:,MOISTURE_IND) all_modes(:,TEMP_IND)];
    flux_file = strcat(flux_save_loc, 'flux', master.setup(loc), '.csv');    
    fid = fopen(char(flux_file),'r+');
    fprintf(fid,'date#,conc1,conc2,conc3,MQ_flux,MM_flux,moisture,temperature');
    fprintf(fid,'\n');
    for line = 1:size(out,1) 
    %Used FGETL to move file pointer a whole line at a time: 
    %See FGETL section below for more information 
        fprintf(fid,'%s,%f,%f,%f,%f,%f,%f,%f', char(dt(line,:)), out(line,1), out(line,2), out(line,3), out(line,4), out(line,5), out(line,6), out(line,7)); 
        fprintf(fid,'\n');
    end; 
    fclose(fid);
  
%}

end


end