clear

wv=[400.0, 425.0, 450.0, 475.0, 500.0, 525.0, 550.0, 575.0, 600.0, 625.0, 650.0, 675.0, 700.0];
wvb=[387.5, 412.5, 437.5, 462.5, 487.5, 512.5, 537.5, 562.5, 587.5, 612.5, 637.5, 662.5, 687.5, 712.5];
nlam=length(wv);
nx=360; ny=160;

fid=fopen('ecco_oasim_edp_387_712_below.bin','r','b');
edp=fread(fid,'float32'); fclose(fid);
edp=reshape(edp, nx, ny,13,12);

fid=fopen('ecco_oasim_esp_387_712_below.bin','r','b');
eds=fread(fid,'float32'); fclose(fid);
eds=reshape(eds, nx, ny,13,12);

ix=20; iy=80; it=4;
pard=squeeze(edp(ix,iy,:,it));
pars=squeeze(eds(ix,iy,:,it));
figure
subplot(2,1,1)
for il=1:nlam,
   plot([wvb(il) wvb(il+1)],[pard(il) pard(il)],'r-','LineWidth',3);
   hold on
   plot([wv(il)],[pard(il)],'ro');
end
subplot(2,1,2)
for il=1:nlam,
   plot([wvb(il) wvb(il+1)],[pars(il) pars(il)],'r-','LineWidth',3);
   hold on
   plot([wv(il)],[pars(il)],'ro');
end


wvn=[400:5:700];
wvbn=[400-2.5:5:700+2.5];
nlamn=length(wvn);


newedp=zeros(nx,ny,nlamn,12);
neweds=zeros(nx,ny,nlamn,12);

for il=1:nlamn,
     clear fi, fi=find(wv==wvn(il));
     if (isempty(fi)==0),
        newedp(:,:,il,:)=edp(:,:,fi,:);
        neweds(:,:,il,:)=eds(:,:,fi,:); 
     else
       clear fi, fi=find(wv<wvn(il));
       fin1=fi(length(fi));
       fin2=fin1+1;
       ax=wvn(il)-wv(fin1);
       bx=wv(fin2)-wvn(il);
       newedp(:,:,il,:)=bx/(ax+bx)*edp(:,:,fin1,:)+ax/(ax+bx)*edp(:,:,fin2,:);
       neweds(:,:,il,:)=bx/(ax+bx)*eds(:,:,fin1,:)+ax/(ax+bx)*eds(:,:,fin2,:);
     end % if
end % for il

pardn=squeeze(newedp(ix,iy,:,it));
parsn=squeeze(neweds(ix,iy,:,it));
subplot(2,1,1)
for il=1:nlamn,
   plot([wvbn(il) wvbn(il+1)],[pardn(il) pardn(il)],'k-','LineWidth',1);
   hold on
%  plot([wv(il)],[pard(il)],'ro');
end
subplot(2,1,2)
for il=1:nlamn,
   plot([wvbn(il) wvbn(il+1)],[parsn(il) parsn(il)],'k-','LineWidth',1);
   hold on
%  plot([wv(il)],[pars(il)],'ro');
end



fid=fopen('ecco_oasim_edp_5nm_below.bin','w','b');
for it=1:12,
 for il=1:nlamn,
  clear tmp, tmp=squeeze(newedp(:,:,il,it));
  fwrite(fid,tmp,'float32');
 end
end
fclose(fid);

fid=fopen('ecco_oasim_esp_5nm_below.bin','w','b');
for it=1:12,
 for il=1:nlamn,
  clear tmp, tmp=squeeze(neweds(:,:,il,it));
  fwrite(fid,tmp,'float32');
 end
end
fclose(fid);


