%Written by Michael.J.S Goh
%Date:24/12/2014
%Custom code to present trend of Continous Mesh CVT
%% Control panel 
ShowRaw=2; %Show raw data  0=no 1=yes 2=show average and raw
EndFit=2; %Type of cleanup for ends, 1=average 2=smooth
fittype='poly4'; %type of fit for shifting section
LB=2; %Start time of gear shift
UB=7; %Stop gear shifting
EL=9; %Simulation End time 
IVel=50; %Input angular velocity  in rpm
filename = 'C:\Users\Michael\Desktop\Cvts\Graphs\Graph 2.4.csv'; %Your data file location
Error=10; % Fluctuation adjustment (not needed, can set to 0. if set to 100 will make almost no adjustment)

graph_title='Gear shift from t=2s to t=7s'; %moved here for convinence

%% Importing the file
delimiter = ',';
startRow = 3;
formatSpec = '%f%f%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'HeaderLines' ,startRow-1, 'ReturnOnError', false);
fclose(fileID);
Time = dataArray{:, 1};
Vel = dataArray{:, 2};
clearvars filename delimiter startRow formatSpec fileID dataArray ans;
Vel=Vel.*60./360./IVel;
Timeo=Time;
Velo=Vel;

%% Flattening the portions where the actuator hub has not moved
Initial=Velo(Timeo(:)<LB);
initial=mean(Initial);
iw=ones(length(Initial),1);
initial=initial.*iw;

Final=Velo(Timeo(:)>UB);
final=mean(Final);
fw=ones(length(Final),1);
final=final.*fw;

%% Processing the velocity while shifting gears
Vel2=Vel(Time(:)>LB);
Time2=Time(Time(:)>LB);
Vel3=Vel2(Time2(:)<UB);
Time3=Time2(Time2(:)<UB);
Vel=Vel3;

error=Error/100;
for i=2:1:length(Vel)-1;
    if Vel(i+1)>(Vel(i-1)+Vel(i)+Vel(i+1))/3*(1+error);
        Vel(i+1)=(Vel(i-1)+Vel(i)+Vel(i+1))/3;
    elseif Vel(i+1)<(Vel(i-1)+Vel(i)+Vel(i+1))/3*(1-error);
        Vel(i+1)=(Vel(i-1)+Vel(i)+Vel(i+1))/3;
    end
end

Time3=[LB.*iw;Time3;UB.*fw];
Vel3=[initial;Vel;Final];

Fitted=fit(Time3,Vel3,fittype);
fitted=Fitted(Time3);

%% Rejoining
if EndFit==1;
    Time4=[Timeo(Timeo(:)<LB);Time3;Timeo(Timeo(:)>UB)];
    Vel4=[initial;fitted;final];
elseif EndFit==2;
    Time4=[0;Time3;EL];
    Vel4=[fitted(1);fitted;fitted(length(fitted))];
end
%% Plotting
if ShowRaw==1;
    plot(Time4,Vel4,'r',Timeo,Velo,'g')
elseif ShowRaw==0;
    plot(Time4,Vel4,'r')
elseif ShowRaw==2;
    plot(Time4,Vel4,'r',Timeo,Velo,'g',Time3,Vel3,'c')
end

xlabel('Time(s)')
ylabel('Output to Input ratio')
title(graph_title);
legend('Cleaned-up','Raw data','Pre-fit')
hold off

