clear all
close all
clc
% PHUA: force rendering __________________________________________
COP_0 = 1;
COP_dz = [0 0.5 1 1.5 2 2.5 3];   % deadzone
A = -2.5; K = 2.5; Q = 0.5;
% Richard's curve: influence of the growth rate B ________________ 
B = 1.5;
COP = linspace(-8,8,101);
%COP = linspace(0,0,101);
F = zeros(length(B),length(COP));

for i=1:1:7
ILog = (COP >= COP_0+COP_dz(i)); % COP > COP0+COP_dz
v = 0.585; M = COP_0+COP_dz(i); 
F(i,ILog) =A+(K-A)./(1+Q*exp(-B*(COP(ILog)-M))).^(1/v);
%I = find(x >= COP_0+COP_dz); I(1), F5(I(1)) % Ajustar F(0)

ILog = (COP < COP_0-COP_dz(i));   % COP < COP0-COP_dz
v = 0.5895; M = COP_0-COP_dz(i); 
F(i,ILog) =A+(K-A)./(1+Q*exp(-B*(COP(ILog)-M))).^(1/v);
%I = find(x < COP_0-COP_dz); I(end), F5(I(end)) % ajustar F(0)



end

 figure
 plot(COP/100 - 0.008,F,'LineWidth',2)
 axis([-0.08 0.08 -2.5 2.5])
 xlabel('COP (mm)','FontSize',14), ylabel('Force (N)','FontSize',14)
 legend('M=0.005','M=0.01','M=0.015','M=0.02','M=0.025','M=0.03','M=0.035')