clear all
%intrduco le variabili
c=299792458; h=6.62607015*10^-34; k=1.380649*10^-23; e=exp(1); %costanti varie ed eventuali
%scrivo le funzioni
flux_calc_full= @(m,temp_star) 6.85*(10^14)*10^(-0.4*m)/temp_star; %calcolo il flusso per tutte freq
areacirc= @(ro, ri) pi*ro^2-pi*ri^2;

dim=linspace(1,0.1,100);
for i=1:length(dim)
    f=@(m) SNR_magn(m,dim(i))-3;
    magn(i)=fzero(f,10);
    f_ten=@(m) SNR_magn(m,dim(i))-10;
    magn_ten(i)=fzero(f_ten,10);
end

res=1.22.*550e-9./dim.*206265;
% flip(dim);
% flip(magn);
% flip(magn_ten);
figure
plot(dim,magn,'r',dim,magn_ten,'g')
colororder({'b','b'})
set(gca, 'YDir','reverse')
title('Maximum Magnitude observable for a given SNR')
xlabel('Diameter of the primary [m]')
ylabel('Apparent Bolometric Magnitude')

%xlim([0 10])
ylim([14 24])
grid on
yyaxis right
hold on
plot(dim,res,'b--')
ylim([0.05 1.5])
ylabel('Diffr. Angular Resolution [arcsec]')

legend('SNR=3','SNR=10','Angular Resolution')
hold off

function SNR = SNR_magn(m,dim)
    exp_time=60; %[s]
    temp=5500; %K
    vis=[400,625]; %[nm]
    %flusso_tot=flux_calc_full(m,temp); %photon flux [m-2 s-1]
    flusso_sens=flux_calc(temp,vis(1),vis(2),m); %[m-2 s-1]
    %fraz_percepita=flusso_sens/flusso_tot;                                %2,73e8 fotoni m-2
%TO DO tieni conto dello scale factor
    areacirc= @(ro, ri) pi*ro^2-pi*ri^2;
    sup_ingresso=areacirc(dim,0.8*dim); %[m2]                                   % @10000k deve essere 0,258
    fluence=flusso_sens*sup_ingresso; %[s-1]
    tr=0.90; %trasmittance mirror alluminio visible
    n_specchi=3; %numero specchi

    QE=0.75; %moltiplico per efficienza quantica
    elec_flux_read=fluence*(tr.^3)*QE; %[electr s-1] %https://scientificimaging.com/knowledge-base/qe-curves-for-cmos-imagers/
    dark_current=3; %[e s-1 p-1] dark current
    read_noise=4; %[e]readount noise

    full_well=20e3; %[e]
%qui dovrei tenere conto della dimensione del disco di airy
    f_number=14;
    centerwl=(vis(2)-vis(1)/2)*10^-9; %[m]
%
    airy=1.22*f_number*centerwl; %[m]
    pixel_size=4*10^-6; %[m]
    pixel_coperti=ceil(airy/pixel_size)^2;
%
    flux_per_pix=elec_flux_read/pixel_coperti;
    CCD_Eq= @(fot)(fot)./sqrt(fot+dark_current*exp_time+read_noise^2); %CCD equasion
    SNR=CCD_Eq(flux_per_pix*exp_time);
end


function f=flux_calc(temp,lamda_s,lamda_l,m) %FLUSSO PER RANGE DI FREQ
    c=299792458; h=6.62607015*10^-34; k=1.380649*10^-23; e=exp(1); %costanti varie ed eventuali
    lamda_s_m=lamda_s*1e-9; lamda_l_m=lamda_l*1e-9; %nm -> m
    df= @(lamda) 1./(lamda.^4.*(e.^(h.*c./(lamda.*k.*temp))-1));
    f_emess=2*pi*c*integral(df,lamda_s_m,lamda_l_m); %fot m-2 superf della stella

    l0=3.828*10^26; %[w] %luminosità del sole
    M0=4.74; %magnitutine bolometrica del sole
    sb=5.670374419*1e-8;
    beta=3.086e16;
    f=f_emess/temp^4*l0/(400*pi*beta^2*sb)*10^(0.4*(M0-m));
end