% view the correlated-k data
close all
clear all

kB = 1.3806504e-23; % boltzmann's constant

vis=0 % in the visible?

cd /u/mvalmd/datagcm/corrk_data/

ls
fstrin=input('Choose dataset to analyse: ')
cd(fstrin)
    
load p.dat
load T.dat
load g.dat

% fid=fopen('Q.dat');
% tmp=fgetl(fid); L_NSPECIES=str2num(tmp)
% for ii=1:L_NSPECIES
%     tmp=fgetl(fid);
%     sname(ii,:)=tmp;
% end
% tmp=fgetl(fid); L_NVAR=str2num(tmp)
% for kk=1:L_NVAR
%     tmp=fgetl(fid);
%     var(kk)=str2num(tmp);
% end
% fclose(fid);
L_NVAR=1

L_NGAUSS  = g(1);
L_NPREF   = p(1);
L_NTREF   = T(1);

Pv=10.^p(2:end);
Tv=T(2:end);
gw=g(2:end)';

clear p T g

ls
fstrin=input('Choose band folder: ')
cd(fstrin)

if(vis)
    fid=fopen('narrowbands_VI.in');
    L_NSPECT=str2num(fgetl(fid));
    for ii=1:L_NSPECT
       bands(ii,:)=str2num(fgetl(fid));
    end
    fclose(fid);
else
    fid=fopen('narrowbands_IR.in');
    L_NSPECT=str2num(fgetl(fid));
    for ii=1:L_NSPECT
       bands(ii,:)=str2num(fgetl(fid));
    end
    fclose(fid);
end

abs4d = zeros(L_NTREF,L_NPREF,L_NVAR,L_NSPECT,L_NGAUSS);
if(vis)
    abs1d = load('corrk_gcm_VI.dat','-ascii');
else
    abs1d = load('corrk_gcm_IR.dat','-ascii');
end
abs4d = reshape(abs1d,[L_NTREF L_NPREF L_NVAR L_NSPECT L_NGAUSS]);


kdistnice=zeros(L_NTREF,L_NPREF);
kdistok  =zeros(L_NTREF,L_NPREF,L_NSPECT);
for iP=1:L_NPREF
    for iT=1:L_NTREF

        kc(:,:)=abs4d(iT,iP,1,:,1:L_NGAUSS);
        kd(1,:)=abs4d(iT,iP,1,1,1:L_NGAUSS);
        %kd(1,:)=abs4d(iT,iP,1,4,1:L_NGAUSS);
        
                   %kc=100^2 * kB * Tv(iT) * kc / ( Pv(iP)  * 100 );
                   %kd=100^2 * kB * Tv(iT) * kd / ( Pv(iP)  * 100 );
                   % convert from k to sigma, SI to cm^2 / molecule
                   
                   kc=Pv(iP)  * 100 * kc / (100^2 * kB * Tv(iT));
                   kd=Pv(iP)  * 100 * kd / (100^2 * kB * Tv(iT));
                   % convert from sigma to k, cm^2 / molecule to SI

        ktemp=sum(kc.*repmat(gw,[L_NSPECT 1]),2)/sum(sum(repmat(gw,[L_NSPECT 1]),2));
    
        for iB=1:L_NSPECT
            kdistok(iT,iP,iB)=ktemp(iB);
        end
        
        kdistnice(iT,iP)=sum(kd.*gw)/sum(gw);
                       
    end
end

figure(1)
subplot(2,1,1)

surf(kdistnice); view(2);
%caxis([0 1e-18])
colorbar horz;  axis  tight

dbands=bands(:,1)-bands(:,2);

figure(2); %hold on
semilogy(mean(bands,2),squeeze(kdistok(end,end,:))','b'); %axis([0 1950 1e-34 1e-18])
cd ..

