// Energy calulation using planck law
h=6.62606957e-34;
k=1.3806488e-23;
c=299792458;
function g=SpectralResponseCCD(L)
    ind0=find(L<400e-9);
    ind1=find(L<600e-9 & L>400e-9);
    ind2=find(L<800e-9 & L>600e-9);
    ind3=find(L<1100e-9 & L>800e-9);
    ind4=find(L>1100e-9 );
    g=zeros(L);
    if length(ind0)>0
        g(ind0)=0;
    end 
    if length(ind1)>0
        g(ind1)=(L(ind1)-400e-9)*0.5/200e-9;
    end 
    if length(ind2)>0
        g(ind2)=0.5;
    end 
    if length(ind3)>0
        g(ind3)=-(L(ind3)-1100e-9)*0.5/300e-9;
    end 
     if length(ind4)>0
        g(ind4)=0;
    end 
endfunction
         
function e=planck(T,lambda)
    e=2*h*c.^2*lambda.^(-5).*(exp(h*c.*(lambda*k*T).^(-1))-1).^(-1);
endfunction

L1=[400:800]*1e-9;// wave band L1 ranging 400 to 800nm
L2=[800:1000]*1e-9;//wave band L2 ranging 800 to 1000nm
e900L1=sum(planck(900,L1).*SpectralResponseCCD(L1))*(L1(2)-L1(1));// energy in L1 with black body at 900K
e900L2=sum(planck(900,L2).*SpectralResponseCCD(L2))*(L2(2)-L2(1));// energy in L2 with black body at 900K
e5800L1=sum(planck(3000,L1).*SpectralResponseCCD(L1))*(L1(2)-L1(1));// energy in L1 with black body at 3000K
e5800L2=sum(planck(3000,L2).*SpectralResponseCCD(L2))*(L2(2)-L2(1));// energy in L1 with black body at 3000K

disp(e900L1/e5800L1) // contrast in L1
disp (e900L2/e5800L2)// contrast in L2 
disp((e900L2/e5800L2)/(e900L1/e5800L1) )// contrast enhancement L2 relative to L1

