demoh=0.06;
demokBT=1.7;
J=1;
FILE_prefix='UMB/FFS_';
lam_pre='lam_array';
tot_pre='total';
FILE_type='.txt';
key=1;
if (abs(round(demoh*100)-demoh*100)<1e-12)
    h_val=sprintf('%4.2f',demoh);
else
    h_val=sprintf('%5.3f',demoh);
end
%h_val=sprintf('0.05');
kBT_val=sprintf(strcat('_',num2str(demokBT),'/'));

% get ncpu
for i=0:3000
    CHECK=strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i),'/Freq.txt');
    key1=exist(CHECK);
    if (key1<1)
        ncpu=i;
        break;
    end
end
Z=zeros(1,ncpu);
Z(1)=1;

% get ncpu
for i=0:4000
    CHECK=strcat(FILE_prefix,h_val,kBT_val,'Kinetic_',num2str(i),'/n_array.txt');
    key1=exist(CHECK);
    kinetic_key=key1;
    if (key1>0)
        FILE1=strcat(FILE_prefix,h_val,kBT_val,'Kinetic_',num2str(i),'/tot_swip.txt');
        temp=load(FILE1);
        Nswip=temp(2);
        Ntot=temp(1);
        
        FILE1=strcat(FILE_prefix,h_val,kBT_val,'Kinetic_',num2str(i),'/time_array.txt');
        time_array=load(FILE1);
        
        FILE1=strcat(FILE_prefix,h_val,kBT_val,'Kinetic_',num2str(i),'/n_array.txt');
        n_array=load(FILE1)/Nswip;
        
        fit_temp=polyfit(time_array,n_array,1);
        kinetic=fit_temp(1)/2;
        
        break;
    end
end


% compute the probability by overlapping window method.

%data=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_arg1',FILE_type));
%Sec_len=data(2)*2+1;
data=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_1/Narray',FILE_type));
Sec_len = length(data);

Narray_TOT=zeros(Sec_len,ncpu);
Freq_TOT=zeros(Sec_len,ncpu);
EXP_TOT=zeros(Sec_len,ncpu);
for i=1:ncpu
    Narraytemp=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-1),'/Narray',FILE_type));
    Freqtemp=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-1),'/Freq',FILE_type));
    Probtemp=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-1),'/Prob',FILE_type));
    
    Narray_TOT(:,i)=Narraytemp;
    Freq_TOT(:,i)=Freqtemp;
    exp_array=Freqtemp./Probtemp;
    exp_array(find(Freqtemp==0))=0;
    EXP_TOT(:,i)=exp_array;
end

for i=2:ncpu
        if (abs(round(demoh*100)-demoh*100)<1e-12)
            h_val=sprintf('%4.2f',demoh);
        else
            h_val=sprintf('%5.3f',demoh);
        end
    Narray1=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(max(i-3,0)),'/Narray',FILE_type));
    Narray2=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-2),'/Narray',FILE_type));
    Narray3=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-1),'/Narray',FILE_type));
    Freq1=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(max(i-3,0)),'/Freq',FILE_type));
    Freq2=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-2),'/Freq',FILE_type));
    Freq3=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-1),'/Freq',FILE_type));
    Prob1=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(max(i-3,0)),'/Prob',FILE_type));
    Prob2=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-2),'/Prob',FILE_type));
    Prob3=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(i-1),'/Prob',FILE_type));
    Z(i)=Z(i-1);
    
   [x,fval] = fminsearch(@sclabs,Z(i),optimset('TolX',Z(i-1)*1e-14),...
       Z(max(1,i-2)),Z(i-1),Narray1,Freq1,Prob1,Narray2,...
            Freq2,Prob2,Narray3,Freq3,Prob3);
    Z(i)=x;
%    i
%    x
%    fval
end

Narraymin=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(0),'/Narray',FILE_type));
Narraymax=load(strcat(FILE_prefix,h_val,kBT_val,'UMB_',num2str(ncpu-1),'/Narray',FILE_type));
n_min=min(Narraymin);
n_max=max(Narraymax);

[aa,bb]=size(Narray_TOT);
p=zeros(1,n_max-n_min+1);
s=0;
for l=n_min:n_max
%    for l=20:20
    ind=find (Narray_TOT == l);
    nume=sum(Freq_TOT(ind));
    deno=sum(  EXP_TOT(ind)'.*...
        sum(Freq_TOT(:,ceil(ind/aa)))./Z(ceil(ind/aa)) );
    s=s+1;
    if(nume==0)
       p(s)=0;
    else
       p(s)=nume/deno;
    end
end

ind2=find(p>10^(-10));
pnorm=sum(p(ind2));
p=p/pnorm;
n_seq=n_min:n_max;

p_new=p(find(p>0));
n_seq_new=n_seq(find(p>0));

figure(701)
semilogy(n_seq_new,p_new)
title('cluster size vs probability');
xlabel('cluster size');
ylabel('probability');

Free_over_kBT=-log(p_new)+log(p_new(1));
figure(702)
plot(n_seq_new,Free_over_kBT,'.-')
title('cluster size vs free energy');
xlabel('cluster size');
ylabel('free energy/kBT');

[aaa,bbb]=min(p_new);
n_c=n_seq_new(bbb);
W_c=demokBT*Free_over_kBT(bbb);

fit_s=ceil(bbb*2/3);
n_seq_new_23=n_seq_new.^(1/2);
[p,S]=...
    polyfit(n_seq_new_23(fit_s:end),...
    demokBT*Free_over_kBT(fit_s:end)+2*demoh*n_seq_new(fit_s:end)...
    -5/4*demokBT*log(n_seq_new(fit_s:end)),1); % Langer correction
Free_fit=p(1)*n_seq_new.^(1/2)-2*demoh*n_seq_new+p(2)...
    +5/4*demokBT*log(n_seq_new);% Langer correction
W_c_new=max(Free_fit);
alpha=p(1);

figure(703)
plot(n_seq_new,demokBT*Free_over_kBT,'b.-',n_seq_new,Free_fit,'k-');

Ntot=10000;

Nc_scalefac=1.00;
hexp=-2;
%demoh=0.40;
%demokBT=2.35;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

J=1;
FILE_prefix='runs/FFS_';
lam_pre='lam_array';
tot_pre='total';
FILE_type='.txt';
%h=[ 0.08 0.13 0.15 0.17 0.20 0.22 0.25:0.05:0.95];
%h=[0.25 0.27 0.30:0.05:0.45 0.50 0.55 0.60 0.70 ];
h=[ demoh ];
%kBT=[ 1.60:0.10:2.30 2.45 2.50:0.10:3.40 ];
kBT = [ demokBT ];
L=0;
DATA=zeros(100,11);
lam=zeros(100,200);
avePb=zeros(100,200);
P0toplus=zeros(100,200);
P0toplus_err=zeros(100,200);

for j=1:length(kBT)
    for i=1:length(h)
        K=J/kBT(j);
        khi=(1-sinh(2*K)^(-4))^(1/8);
        sig_parr=kBT(j)*(2*K+log(tanh(K)));
        sig_diag=kBT(j)*(sqrt(2)*log(sinh(2*K)));
        sig_eff=(sig_parr+sig_diag)/2/sqrt(khi);
        M=(1-1/sinh(2/(kBT(j)))^4 )^(1/8);
        Nc=pi*(sig_eff/(h(i)*2))^2/M;
        Nc_sol=Nc;
        if (demokBT==kBT(j))
            sig_se=sig_eff;
            x_se=M;
        end
        Wc=pi*(sig_eff^2)/(h(i)*2);         
        Wc_sol=Wc;
        if (abs(round(demoh*100)-demoh*100)<1e-12)
            h_val=sprintf('%4.2f',demoh);
        else
            h_val=sprintf('%5.3f',demoh);
        end
        kBT_val=sprintf('_%3.1f/',kBT(j));
        
        %computing predictions

        CHECK=strcat(FILE_prefix,h_val,kBT_val,'autoend.txt');
        key1=exist(CHECK);
        P0to_temp=zeros(1,200);
        P0to_temp_err=zeros(1,200);
        if (key1>0) 
            for k=0:100
                no=sprintf('%d',k);
                FILE=strcat(FILE_prefix,h_val,kBT_val,tot_pre,no,FILE_type);
                key=exist(FILE);
                if (k==0 && key ==0)
                    break;
                elseif (k==0 && key >0)
                    L=L+1;
                elseif (k>0 && key==0)
                    break;
                end
                
                data=load(FILE);
                if (k==0)
                    initN=data(2);
                    Init_Flux=data(2)/data(1)/Ntot;
                elseif (k==1)
                    Pcon=data(2)/data(1);

                    err_cur=sqrt(Pcon*(1-Pcon)/data(1));
                    err_cum=err_cur;
                    Pcum=Pcon;
                else
                    Pcon=Pcon*data(2)/data(1);
                    Pcur=data(2)/data(1);

                    err_cur=sqrt(Pcur*(1-Pcur)/data(1));
                    err_cum=Pcon*sqrt(err_cur^2/Pcur^2+err_cum^2/Pcum^2);
                    Pcum=Pcon;
                end
                
                if (k>=1)
                    P0to_temp(k)=Pcon;
                    P0to_temp_err(k)=err_cum;
                end
                    
                prev_k=k;
            end
            TOT=prev_k;
            k_AB=Init_Flux*Pcon;
            AB_err=err_cum/Pcon;
            
            LAM=strcat(FILE_prefix,h_val,kBT_val,'lam_array.txt');
            
            fp=fopen(LAM);
            fgets(fp);fgets(fp);
            lam_temp=zeros(1,200);
            for l=1:200               
                lam_temp(l)=str2num(fgets(fp));
            end
            fclose(fp);
            
            lam(L,:)=lam_temp;
            P0toplus(L,:)=P0to_temp;
            P0toplus_err(L,:)=P0to_temp_err;     
            for l=TOT:-1:1
                no=sprintf('%d',l);
                FILE=strcat(FILE_prefix,h_val,kBT_val,tot_pre,no,FILE_type);
                data=load(FILE);
                cond_P=data(2)/data(1);
                if (l<TOT)
                    avePb(L,l)=cond_P*avePb(L,l+1);
                else
                    avePb(L,l)=cond_P;
                end
            end

            Nc_pb=0;
%            TOT=Final_data(L,8);

            [a,b]=min(abs(avePb(L,1:TOT)-0.3));
            [c,d]=min(abs(avePb(L,1:TOT)-0.7));            
            if (b < d )
                fito=max(d-b-2,1);
                if (fito>5)
                    fito=4;
                end
                disp(sprintf('kBT=%f, h=%f d-b=%d fito=%d initN=%d',...
                    kBT(j),h(i),d-b,fito,initN));
                coef=polyfit(lam(L,b:d),avePb(L,b:d),fito);
                xaxis=lam(L,b):lam(L,d);
                yaxis=polyval(coef,xaxis);
                [a,cind]=min(abs(yaxis-0.5));
                Nc_pb=xaxis(cind);
            else
                disp(sprintf('kBT=%f h=%f, dont believe even pb!\n',...
                    kBT(j),h(i)));
            end
            
            if (avePb(L,TOT) < 0.5 )
                Nc_pb=0;
            end
            
            DATA(L,:)=[ kBT(j) h(i) Init_Flux Pcon k_AB Wc Nc TOT ...
                Nc_pb avePb(L,TOT) AB_err];
        end        
    end
end

Final_data=zeros(L,5);
Final_data=DATA(1:L,:);
Final_data;

%h dependence
ratiosum=0;
ratiosum_l_Tr=0;
ratiosum_b_Tr=0;
ratio_array=zeros(length(kBT),1);

% Draw avePb plot
plot_no=i+j+1;

index=find(Final_data(:,1)==demokBT & Final_data(:,2)==demoh);
TOT=Final_data(index,8);
[a,b]=min(abs(avePb(index,1:TOT)-0.2));
[c,d]=min(abs(avePb(index,1:TOT)-0.8));
if (b < d )
    fito=max(d-b-1,1);
    if (fito>5)
        fito=4;
    end
coef=polyfit(lam(index,b:d),avePb(index,b:d),fito);
xaxis=lam(index,b):lam(index,d);
yaxis=polyval(coef,xaxis);
[Nc_pb,cind]=min(abs(yaxis-0.5));
Nc=xaxis(cind)
figure(plot_no)
plot(lam(index,1:TOT),avePb(index,1:TOT),'o',xaxis,yaxis);
else
    disp(sprintf('kBT=%f h=%f, do this again!',kBT,h));
end
title(sprintf('cluster size vs free energy',demokBT,demoh))
xlabel('cluster size');
ylabel('free energy/kBT');


% Draw P0toplus plot
plot_no=plot_no+1;

index=find(Final_data(:,1)==demokBT & Final_data(:,2)==demoh);
TOT=Final_data(index,8);
figure0=figure(plot_no);
axes('Parent',figure0,'YScale','log','YMinorTick','on');
box('on');
hold('all');
errorbar(lam(index,2:TOT),P0toplus(index,1:TOT-1),...
    P0toplus_err(index,1:TOT-1),'o-');
title(sprintf('Conditional grow probability at kBT = %f & h = %f',demokBT,demoh))
xlabel('Size of largest nucleus')
ylabel('P(\lambda_{i+1}|\lambda_i)')
hold off;

% compare FFS and UMB
P_ref=p_new(lam(index,1)+1);
P_FFT=P_ref*P0toplus(index,1:TOT-1);
figure(801)
semilogy(lam(index,2:TOT),P_FFT(index,1:TOT-1),'o',n_seq_new,p_new,'-');

disp(sprintf('\n'));
disp(sprintf('h=%f kBT=%f',demoh,demokBT));
disp(sprintf('n_c=%d',n_c));
disp(sprintf('W_c=%f',W_c));
disp(sprintf('h=%f, W_c/n_c=%f',demoh,demokBT*Free_over_kBT(bbb)/n_c));
Zeldo=sqrt(W_c/(4*pi*demokBT*n_c^2)); % BD CNT

n_ccc=(( (alpha/2)+sqrt( alpha^2/4+8*h*1.25*kBT ))/(4*h))^2;

%xiaoha: two undefined variables
n_c_fit = 1; k_AB = 1;

lamL=0.25*alpha*n_c_fit^(-3/2)+1.25*kBT/(n_c_fit^2);% with Langer Correction
Zeldo=(lamL/(2*pi*demokBT))^(1/2);

disp(sprintf('Zeldovich Factor=%f',Zeldo));
[ac,bc]=max(Free_fit);
n_c_fit=n_seq_new(bc);
disp(sprintf('sigma=%f n_c=%f Wc=%f W_c/n_c=%f from Fit\n',p(1),...
    n_c_fit,W_c_new,W_c_new/n_c_fit));

if (kinetic_key>0)

figure(805)
plot(time_array,n_array);
title('timestep vs (\Delta n)^2');
xlabel('timestep');
ylabel('(\Delta n)^2');

I_BD=Zeldo*kinetic*exp(-W_c/demokBT);
disp(sprintf('kinetic = %10.5e',kinetic));
disp(sprintf('K_BD = %10.5e',I_BD));
disp(sprintf('K_AB = %10.5e',k_AB));
disp(sprintf('K_BD/k_AB = %10.5f',I_BD/k_AB));
end
disp(sprintf('Langer_prefac = %10.5f',k_AB/(exp(-W_c/demokBT))));
disp(sprintf('\nFrom 2D Solution \nNc=%10.5f Wc=%10.5f coeff=%10.5f',...
    Nc_sol,Wc_sol,2*sqrt(pi)*sig_eff/sqrt(M)));

W_c0=n_c_fit*2*h;
W_c1=W_c0+1.25*kBT*log(n_c_fit);
prefact=k_AB/(  exp(-W_c0/kBT) );
dWc=W_c-W_c0;
disp(sprintf('W_c0=%10.5f, W_c1=%10.5f dWc=%10.5f exp(-dWc/kBT)=%f',...
    W_c0,W_c1,dWc,exp(-dWc/kBT)));
disp(sprintf('prefact=%10.5e',prefact));
