%input variables
force = input('Applied Force: ');
day = input('Number of Days: ');

a = 1;
b = 2;
counter = 1;

%determine stress
area = 100; %area in mm^2
stress = force / area;

%given constants!!
kD = 1.85 * 10^5; %damage rate coefficient in mm/mm^2
RL = 3000; %cycles of loading applied per day
A0 = 0.02835; %area of a BMU in mm^2
Fs = 5; %repair specificity factor
famax = 0.5; %BMUs/mm^2/day
kr = -1.6;
%NF = 0.4288; BMUs/mm^2
QR = 0.00113; %mm^2
QF = 0.000423; %mm^2
PHC= 0.0433;

Porosity(1) = 0.0500; %initial porosity

%Activation Frequency, NR, NF
for a = 1:95

    ActivationFrequency(a) = 0.00670; %initial activation frequency in BMUs/mm^2/day

   if a <= 25
    NR(a) = sum(ActivationFrequency(1:a));

   elseif   25 < a && a < 30
        NR(a)= NR(25);
        NF(a) = 0;

   elseif 30 <= a && a <= 94
        NF(a) = sum(ActivationFrequency(1:a));
        NR(a) = NR(25);
    end

end

%Porosity
for b = 2:95
    c= Porosity(b - 1);
    d = ((c) + (((1 - PHC) * QR * NR) - (QF * NF)));
    Porosity(b)= d;

    E = 23400 * (1 - Porosity(b - 1))^5.74;  %E in MPA

    strain = stress / E;

    Dformed = kD * RL * strain^4; %damage formed in mm/mm^2

    Dremoved = Damage(b - 1) * ActivationFrequency(b - 1) * A0 * Fs;

    Damage(b) = Damage(b - 1) + Dformed - Dremoved; %damage in mm/mm^2
end


%for counter= 2:95

   % Porosity(counter) = Porosity(counter - 1) + (((1 -PHC) * QR * NR) - (QF * NF));
   % E = 23400 * (1 - Porosity(counter - 1))^5.74;  %E in MPA

   % strain = stress / E;

  %  Dformed = kD * RL * strain^4; %damage formed in mm/mm^2

  %  Dremoved = Damage(counter - 1) * ActivationFrequency(counter - 1) * A0 * Fs;

   % Damage(counter) = Damage(counter - 1) + Dformed - Dremoved; %damage in mm/mm^2

   % ActivationFrequency(counter) = fa0;
   % NR = sum(ActivationFrequency(25:(counter))); %NRprevious + fa;
   % NF = sum(ActivationFrequency(94:30));
   % Porosity(counter) = Porosity(counter - 1) + (((1 -PHC) * QR * NR) - (QF * NF));
    %determine bone density Density(counter) = 2 * (1 - P(counter)); %density in g/cm^3

  %  counter = counter + 1;
%end

for counter= 96:(day + 95)

    %determine damage (D)
    E = 23400 * (1 - Porosity(counter - 1))^5.74;  %E in MPA

    strain = stress / E;

    Dformed = kD * RL * strain^4; %damage formed in mm/mm^2

    Dremoved = Damage(counter - 1) * ActivationFrequency(counter - 1) * A0 * Fs;

    Damage(counter) = Damage(counter - 1) + Dformed - Dremoved; %damage in mm/mm^2

    %determine activation frequency
    ActivationFrequency(counter) = (fa0 * famax)/(fa0 + (famax - fa0) * exp((kr * famax * ((Damage(counter) - D0)/D0))));

    %determine bone density Density(counter) = 2 * (1 - P(counter)); %density in g/cm^3
    NR = sum(ActivationFrequency(25:(counter))); %NRprevious + fa;
    NF = sum(ActivationFrequency(94:30));
    Porosity(counter) = Porosity(counter - 1) + (((1 -PHC) * QR * NR) - (QF * NF)); %Pprevious + ((QR * NR) - (QF * NF));

    counter = counter + 1;

end

%determine porosity (P)




for i = 1:95

    Counter2(i) = i;

end

figure(1)
plot(Counter2,Damage)