clc
clear all
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Be Name Khoda ! &&&&
 % Mesal 3-3   
 %%%%%%%%%%%%%%%%%%%%%%%
 m=1000;
 
 MW0= 72 %pcf
 sigma_MW=.05.*MW0;
MW=normrnd(MW0,sigma_MW,[m 1]);
 
 
 
Ps0 = 0.1 % psi / ft % pheshar sayal
 sigma_Ps=.05.*Ps0;
Ps=normrnd(Ps0,sigma_Ps,[m 1]);


Pr0 = 8000 % psi Pheshar_Makhzan 
 sigma_Pr=.05.*Pr0;
Pr=normrnd(Pr0,sigma_Pr,[m 1]);

Lc0 = 3600 %Meter toole loole casing
 sigma_Lc=.05.*Lc0;
Lc=normrnd(Lc0,sigma_Lc,[m 1]);

M=2250
 N=1650
 P=0
 
 noghte = [M N P]
 BSFT = [6360 5760 5760 6360]
 disp('Noghte =M ya Noghte =N Ya Noghte P ! ')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 disp(' Mohasebat Dar Noghte')
     disp('M')
     disp('==================================================')
 P1= Pr - ( (M).* 3.281.* Ps )
 P0=(MW./144).*3.281.*M
 BL= P1 - P0
 BSF1 = BSFT(1) ./ BL
 
 TR=linspace(1.2,3.1,50);

ss=zeros(size(TR));
pp=zeros(size(TR));

for i=1:length(TR)
  ss(i)=sum(BSF1>=TR(i));  
  pp(i)=ss(i)./m*100;
end

figure(1)
plot(TR,pp);
xlabel('BSF1');
ylabel('Exceedance Probability');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 disp(' Mohasebat Dar Noghte')
     disp('N')
     disp('==================================================')
 P1= Pr - ( (N).* 3.281.* Ps )
 P0=(MW./144).*3.281.*N
 BL= P1 - P0
 BSF2 = BSFT(2) ./ BL
 
 TR=linspace(0.9,1.7,50);

ss=zeros(size(TR));
pp=zeros(size(TR));

for i=1:length(TR)
  ss(i)=sum(BSF2>=TR(i));  
  pp(i)=ss(i)./m*100;
end

figure(2)
plot(TR,pp);
xlabel('BSF2');
ylabel('Exceedance Probability'); 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 disp(' Mohasebat Dar Noghte')
     disp('P')
     disp('==================================================')
 P1= Pr - ( (P).* 3.281.* Ps )
 P0=(MW./144).*3.281.*P
 BL= P1 - P0
 BSF3 = BSFT(3) ./ BL
 
 TR=linspace(0.6,0.9,50);

ss=zeros(size(TR));
pp=zeros(size(TR));

for i=1:length(TR)
  ss(i)=sum(BSF3>=TR(i));  
  pp(i)=ss(i)./m*100;
end

figure(3)
plot(TR,pp);
xlabel('BSF3');
ylabel('Exceedance Probability');
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 Mg=0.5 % mud gradient
 Gg= 0.1 %gas gradient
 BSf=1.1
 BP=5760 %psi
 WHP=Pr-Lc.*3.281.*Gg % well head pressure
 L=(WHP./(Mg-Gg))-(BSFT(2)./(BSf.*(Mg-Gg))) % tool loole bar hasbe foot
 
Lave_ft=mean(L);
Lave_m=mean(L)/3.281;
figure(4)
hist(L/3.281,50)

BSF4 = BSFT(4)./WHP 
TR=linspace(0.7,1.2,50);

ss=zeros(size(TR));
pp=zeros(size(TR));

for i=1:length(TR)
  ss(i)=sum(BSF4>=TR(i));  
  pp(i)=ss(i)./m*100;
end

figure(5)
plot(TR,pp);
xlabel('BSF4');
ylabel('Exceedance Probability');



%  if BSF>=1 
%      disp(' Monaseb ')
%  else
%      disp(' Na Monaseb')
%  end
%  
%  for i=1:length(noghte)
%      disp(' Mohasebat Dar Noghte')
%      disp(noghte(i))
%      disp('==================================================')
%  P1= Pheshar_Makhzan - ( (tool-noghte(i)) * 3.281* Ps )
%  P0=(MW0/144)*3.281*noghte(i)
%  BL= P1 - P0
%  BSF = BSFT(i) / BL
%  if BSF>=1 
%      disp(' Monaseb ')
%  else
%      disp(' Na Monaseb')
%  end
%  end