clear all; close all;
CurrentFolderPath = pwd;
addpath([CurrentFolderPath,'\common\']);

%% User Inputs


MissionDir = 'C:\jvdata\2018test\';
FlxFileName = '160806_ATOM1_z01_firstdecent_FLX_280-640nm.tab';  

UseTAP = 1;
TapFileName = '160806_jv.tap';

% MissionDir = 'C:\jvdata\2016ATOm1\';
% FlxFileName = '160806_ATOM1_z01_firstdecent_FLX_280-640nm.tab';  
% 
% UseTAP = 1;
% TapFileName = '160806_jv.tap';

StartWvlCutoff = 298;
EndWvlCutoff = 640;

RxnListFileName = 'RXN_List_MatlabTUV53_all.txt'
%% Set File paths
FlxFilePrefix = strtok(FlxFileName,'.');
OutputFileName = [FlxFilePrefix,'_MATLAB_TUV53.jv'];

FlxFilePath = [MissionDir,'flx\',FlxFileName];

JvOutPath = [MissionDir,'jv\',OutputFileName];

XSQY_ProgramDir = [CurrentFolderPath,'\Fortran\'];
XSQY_ProgramName =  'XSQY_TUV53';


XSQY_ProgramPath = [XSQY_ProgramDir,XSQY_ProgramName,' MatlabDriver.txt'];
XSQY_OutputPath = [XSQY_ProgramDir,'\XSQYtest.txt'];

RxnListFile = [pwd,'\common\',RxnListFileName];

%% Read Flux
%Wvl = [300 325 350 375 400 425 450 475 500 525 550 575 600 625 650];
[FluxOut,TimeFlx,Wvl,Header] = ReadFlxTab(FlxFilePath);
nTimeFlx = length(TimeFlx.SAM);
nWvl = length(Wvl);

% calculate flux bin widths from Wvl (center wavelength of pixel)
[discard,iStartWvl] = min(abs(Wvl - StartWvlCutoff));
[discard,iEndWvl] = min(abs(Wvl - EndWvlCutoff));

StartWvl = Wvl(iStartWvl);
EndWvl = Wvl(iEndWvl);

Wvl = Wvl(iStartWvl:iEndWvl);
FluxOut = FluxOut(:,iStartWvl:iEndWvl);
nWvl = length(Wvl);


BinWidth = zeros(nWvl,1);

BinWidth(1) = Wvl(2) - Wvl(1);
for iWvl = 2:nWvl-1
   BinWidth(iWvl) = .5 * (Wvl(iWvl+1) - Wvl(iWvl-1));
end
BinWidth(nWvl) = Wvl(end) - Wvl(end-1);

% Run dummy iteration to get nRXN's

TempC1 = 0;
PressureMB1 = 990;

nRxn = WriteXSQYDrvrFile(Wvl,TempC1,PressureMB1,XSQY_ProgramDir,RxnListFile);
RunFortranCode;  %Fortran in separate script to get 'profiler' process TimeFlx
%ReadXsqyFile;
[Wvl,XSQY,jlabel,RxnNumber] = ReadXsqyTUV53(XSQY_OutputPath,nRxn);

% Temperature and Pressure for JV calculations
if UseTAP == 1  
    TapFilePath = [MissionDir,'tap\',TapFileName];
    [TimeTAP,PressMB,TempC] = ReadJvTap(TapFilePath);   %Reads 3 column file with TimeFlx(Seconds after midnight, tempC, pressMB);
    nTimeTAP = length(TimeTAP);

    %remove bad data points 
    for i = 1:nTimeTAP
        if TempC(i) < -150 || TempC(i) > 150
            TempC(i) = -151;
        end
        if PressMB(i) > 1100 || PressMB(i)< 1;
            PressMB(i) = -9;
        end
    end
    
    iGoodT = find(TempC~= -151);
    iGoodP = find(PressMB~= -9);
    
    TempC_interp = interp1(TimeTAP(iGoodT),TempC(iGoodT),TimeFlx.SAM,'linear',-9);  %set to -9999 if outside time range
    PressMB_interp = interp1(TimeTAP(iGoodP),PressMB(iGoodP),TimeFlx.SAM,'linear',-9);
 
    f1 = figure(1) 
    set(f1,'Position',[100 100 800 800])
    subplot(2,1,1)
    plot(TimeTAP/3600,TempC,'or');
    hold on
    plot(TimeFlx.SAM/3600,TempC_interp,'-b','Linewidth',2)
    hold off
    ylabel('Temperature (C)')
    grid on
    legend('TAPFile','Interp')
    
    subplot(2,1,2)
    plot(TimeTAP/3600,PressMB,'or');
    hold on
    plot(TimeFlx.SAM/3600,PressMB_interp,'-b','Linewidth',2)
    hold off
    ylabel('Pressure (mb)')
    xlabel('Time UTC(hrs)')
    grid on
    legend('TAPFile','Interp')
    
else  
    % set to standard or user defined pressure/temperature
    TempC = ones(nTimeFlx,1)*298;
    PressMB = ones(nTimeFlx,1)* 1013;
end


% Open Output file File 
%Write headers on output
disp(['Writing Output: ',JvOutPath]);
fidOut = fopen(JvOutPath,'w');
fprintf(fidOut,'6\n');

fprintf(fidOut,'JV processing version: MATLAB_TUV53_180207\n');
fprintf(fidOut,'Filler\n');
fprintf(fidOut,'Filler\n');
fprintf(fidOut,'Filler\n');
fprintf(fidOut, '%s\n',['FluxInputFile: ',FlxFilePath])
fprintf(fidOut,'Filler\n');
fprintf(fidOut,'%s%0.2f%s%0.2f\n','JV Start Wvl: ',StartWvl,', JV End Wvl: ',EndWvl);

    
fprintf(fidOut,'%i\n',nRxn);


for iRxn = 1:nRxn
   fprintf(fidOut,'\t%s',jlabel(iRxn).txt);
end
fprintf(fidOut,'\n');

tic    
JvData = zeros(nTimeFlx,nRxn);
for iTime = 1:nTimeFlx
%for iTime = 1:15
%     if mod(iTime,100) ==0
%         disp(['TimeFlx ',num2str(iTime),'/',num2str(nTimeFlx)]);

%     end
    
    disp(['TimeFlx ',num2str(iTime),'/',num2str(nTimeFlx)]);

    
    if TempC_interp(iTime) == -151 || PressMB_interp(iTime) == -9
        JvData(iTime,:) = -9999; 
    else
        nRxn = WriteXSQYDrvrFile(Wvl,TempC_interp(iTime),PressMB_interp(iTime),XSQY_ProgramDir,RxnListFile);
        RunFortranCode;  %Fortran in separate script to get 'profiler' time
        %ReadXsqyFile;
        [Wvl,XSQY,jlabel,RxnNumber] = ReadXsqyTUV53(XSQY_OutputPath,nRxn);
        JvData(iTime,:) = CalcJV(Wvl,BinWidth,FluxOut(iTime,:),XSQY);
    end
    fprintf(fidOut,'%.1f',TimeFlx.SAM(iTime));
    fprintf(fidOut,'\t%.3g',JvData(iTime,:));   %for each Rxn
    fprintf(fidOut,'\n');
end
toc


fclose(fidOut);
disp(['File Written: ',JvOutPath]);
%%
iJO3 = 1;
iJNO2 = 3;

f2 = figure(2)
set(f2,'Position',[700 100 600 800])
subplot(2,1,1)
plot(TimeFlx.SAM,JvData(:,iJNO2),'ob')
grid on
ylabel('jNO2')

subplot(2,1,2)
plot(TimeFlx.SAM,JvData(:,iJO3),'ob')
grid on
ylabel('jO3')