clc
close all
clear all

%%%%%%%%%%%%%%%% Soal 
disp( ' % N-80 , 72# ' )
m=1000;
ra = 10; % Allowable risk
tl=2e3;
sigma_tl=.05.*tl;
tool_loole= normrnd(tl,sigma_tl,[m 1]);

figure(1)
hist(tool_loole,50)
xlabel('Casing length')
ylabel('frequency')

CR=2670;     % collaps Resistance psi

MW0=75;       %(pam)
sigma_MW=.05.*MW0;
MW=normrnd(MW0,sigma_MW,[m 1]);

figure(2)
hist(MW,50)
xlabel('Mw')
ylabel('frequency')

W=72;         % vazne vahed tool (lb/ft)
A=20.7677;    % sathe maghtae loole bala(loole zaaiftar)
K=2.7296;
Max_Hole_P = tool_loole.* 3.281.* MW./ 144; % (psi)
Dx1= ( CR.* 144)./(3.281.*MW); % omgh barasi (m)
BF = 1 - MW./489.5;                   % Buoyancy factor
Wb = (2000 - Dx1).*3.281.*W.*BF;    % buoyancy weight

Sa=(Wb./A) - (K.*( Dx1.*3.281.*MW./144 ) ); % Axial stress
P_Sa=sum(Sa<0)/m*100;

disp(' ------------------------------------------ ')
%%%%%%%%%%%%%%%%%%%%%%%%
% J-55 , 68#
CR=1950;
A=19.4451;
K=2.8676;

ym0=55000; % minimum yield stength
sigma_ym=.05*ym0;
ym=normrnd(ym0,sigma_ym,[m 1]);

figure(3)
hist(ym,50)
xlabel('ym')
ylabel('frequency')

TCR=1950;  % Tabulate Collapse Resistance
Dx2=(TCR.* 144)./( 3.281.*MW); % omgh barasi badi
Wb1=Wb + ( (Dx1-Dx2).*3.281.*W.*BF);

Sa1=(Wb1./A) - (K.*( Dx2.*3.281.*MW./144 ) ); % axial Stress
L2=sum(Sa1<0);
P_Sa1=L2/m*100

X1=Sa1./ ym; % zarib etminan
CCF1 =  abs( 0.5.*( X1-sqrt(4-3.*X1.*X1 ))); % Collapse Correction Factor
CCR1 = TCR.* CCF1 ; % Corrected Collapse Resistance

CSF1= CCR1./ (Dx2.*3.281.*MW./144 ); % collapse Saftey Factor
% TR=[.7 .75 .8 .85 .9 .95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4];
TR=linspace(.9,1,50);

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

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

figure(4)
plot(TR,pp);
xlabel('CSF1');
ylabel('Exceedance Probability');

%CSF_S1 = interp1(pp,TR,ra)

CSF_S1 = input('If allowable risk = 10, CSF = ');



if(CSF_S1<0.999)
    disp(' CSF1<1 ==> Pas Bayad Omgh TasHih Shavad ! ');
    Dx3= ( CCR1.* 144)./ (3.281.*MW);
end

if(CSF_S1>0.999)
    disp(' Ghabel Qabool !  ');
    return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Wb2= 87639 + ( (Dx1-Dx3).*3.281.*W.*BF);
Sa2=(Wb2./A) - (K.*( Dx3.*3.281.*MW/144 ) );
L3=sum(Sa2<0);
P_Sa2=L3/m*100;


X2=Sa2./ ym;
CCF2 =  abs( 0.5.*( X2-sqrt(4-3.*X2.*X2 ))); % Collapse Correction Factor
CCR2 = TCR.*CCF2;      % Corrected Collapse Resistance
CSF2= CCR2 ./ (Dx3.*3.281.*MW./144 ); % collapse Saftey Factor

TR=linspace(.9,1.1,50); 
ss=zeros(size(TR));
pp=zeros(size(TR));

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

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

%CSF_S2 = interp1(pp,TR,ra)
CSF_S2 = input('If allowable risk = 10, CSF = ');

if(CSF_S2<0.999)
    disp(' CSF2<1 ==> Pas Bayad Omgh TasHih Shavad ! ');
    Dx4= ( CCR2.* 144)./ (3.281.*MW);
end

if(CSF_S2>0.999)
    disp(' Ghabel Qabool !  ');
    return;
end

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

Wb3= 87639 + ( (Dx1-Dx4).*3.281.*W.*BF);
Sa3=(Wb3./A) - (K.*( Dx4.*3.281.*MW/144 ) );
L4=sum(Sa3<0);
P_Sa3=L4/m*100;


X3=Sa3./ ym;
CCF3 =  abs( 0.5.* ( X3-sqrt(4-3.*X3.*X3 )));
CCR3 = TCR.* CCF2;
CSF3= CCR3./ (Dx4.*3.281.*MW./144 );

TR=linspace(.9,1.2,50); 
ss=zeros(size(TR));
pp=zeros(size(TR));

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

figure(6)
plot(TR,pp);
xlabel('CSF1');
ylabel('Exceedance Probability');

%CSF_S3 = interp1(pp,TR,ra)
CSF_S3 = input('If allowable risk = 10, CSF = ')

if(CSF_S3<0.999)
    disp(' CSF<1 ==> Pas Bayad Omgh TasHih Shavad ! ')
    Dx5= ( CCR3.* 144)./ (3.281.*MW);
end

if(CSF_S3>0.999)
    disp(' Ghabel Qabool !  ')
    return;
end