clear
A=100;
KD=185000;
RL=3000;
AO=.02835;
FS=5;
QR=.00113;
QF=.000423;
NF=.4288;
FaMax=0.5;
KR=-1.6;
delta_t=1;
DO=.0366;
FAO=.00670;
PHC=0.0433;

D(1:95)=.0366;
FA(1:95)=.00670;
NR(1:95)=.1675;
NF(1:95)=.4288;
P(1:95)=.0500;
density(1:95)=1.90;

a=[FA(1)];

t=input('Number of days for analysis =')+95;
F=input('Force in Newtons = ');

for i=95:t

    E=23400*(1-P(i))^5.74;
    stress = F/A;
    strain = stress/E;

    Dformed=KD*RL*(strain)^4;
    Dremoved=D(i)*FA(i)*AO*FS;

    D(i+1)=D(i)+Dformed-Dremoved;

    FA(i+1)=(FAO*FaMax)/(FAO+(FaMax-FAO)*exp(KR*FaMax*((D(i+1)-DO)/DO)));

    NR(i+1)=sum(FA(i-24:i+1))*delta_t;
    NF(i+1)=sum(FA(i-93:i-29))*delta_t;

    P(i+1)=P(i)+((1-PHC)*QR*NR(i+1)-QF*NF(i+1));

    density(i+1)=2*(1-P(i+1));

    a=[a;FA(i+1)];

end

figure(1)
plot(D)
xlabel('time(days)')
ylabel('Damage (mm/mm^2)')
title('Bone Growth')
format('long');

figure(2)
plot(density)
xlabel('time *days)')
ylabel('density(g/cm^3)')
title('Density of Bone')

disp ('D=');
disp(D(i+1));