function[input]=SaturationCorrectionVector(output,binwid,binnum,shotsnum);

% PMT Saturation Correction using Newtonian method
% For Maui Na lidar, discriminator dead time is 10 ns
% and "PMT + discriminator" maximum count rate is 55 MHz. 

global pi c Me Qe E0 Oss kB NA AMU Radius freqi strength
global FreqNaD2a FreqNaCrossOver FreqAOshift Lambda_center SigmaRay
global Deadtime_Disc No_of_Freqs PDA_Freq_Offset SigmaL MAX_Countrate

taoD=Deadtime_Disc*1e6;                % Discriminator dead time (us)
taoP=(1/MAX_Countrate-Deadtime_Disc)/exp(1)*1e6;     % PMT dead time (us)
tbin=0.32*binwid/.048;						% time for each bin
eps=0.005;										% count rate precision
x0(1:binnum)=1.5;											% initial value for iteration
y=output/shotsnum/tbin;								% output photon cnts ->output count rate

x=x0-(x0.*exp(-taoP*x0)./(1+x0*taoD.*exp(-taoP*x0))-y)./exp(-taoP*x0)./(1-taoP*x0).*(1+taoD*x0.*exp(-taoP*x0)).^2;
while (~isempty(find(abs(x.*exp(-taoP*x)./(1+x*taoD.*exp(-taoP*x))-y)>eps)) | ~isempty(find(abs(x-x0)>eps)))
x0=x;
x=x0-(x0.*exp(-taoP*x0)./(1+x0*taoD.*exp(-taoP*x0))-y)./exp(-taoP*x0)./(1-taoP*x0).*(1+taoD*x0.*exp(-taoP*x0)).^2;
end
input=x*shotsnum*tbin;							%  input count rate -> input photon cnts	

% Newtonian method is used in the saturation correction ----
% The nonlinear equation is f(x)=x*exp(-taoP*x)/(1+x*taoD*exp(-taoP*x))-y=0
% The approximate solution of the equation is x=x0-f(x0)/f'(x0), 
% where x0 is the initial value of the solution, and this solution will be
% continued x1=x0-f(x0)/f'(x0), x2=x1-f(x1)/f'(x1), ...
% until x*exp(-taoP*x)/(1+x*taoD*exp(-taoP*x))=y and x=x0 within the error limit eps.