% calculate the hysteresis curve of the Bastelmotor
% from current measurements in the field coil and
% voltage measurements of the rotor coil

% 2024-07-10 ch

data=csvread("data.csv");
start=3;
t=data(start:end,1);
iField= data(start:end,2); # current sensor voltage
vField= data(start:end,3); # voltage applied to field coil
vRotor= data(start:end,4); # voltage measured on rotor coil
vAmp= data(start:end,5); % voltage measured at poser amplifier output

subplot(3,1,1)

plot(t,iField,"y")
hold on
plot(t,vField,"m")
plot(t,vRotor,"c")
plot(t,vAmp,"g")
title("raw signals")
legend("iField","vField","vRotor","vAmp");
grid on
hold off

# calculate physical values
# ACS712x05B sensor voltage to current relation
acs712_mV_per_A=185;
c=acs712_mV_per_A/1000;
# current calculated from sensor voltage
# negative current to correct hysteresis plot
iField_A=-iField/c;
Rshunt=5.1
# current calculated from difference voltage over the shunt resistor
iFieldR_A=-(vAmp-vField)/Rshunt;
xlabel("time [s]")
ylabel("u [V]")

subplot(3,1,2)
plot(t,iField_A);
hold on 
plot(t,iFieldR_A)
grid on
title("current from sensor and current from shunt");
legend("sensor","Rshunt");
hold off
xlabel("time [s]");
ylabel("i [A]");

subplot(3,1,3)
plot(t,vRotor);
grid on
title("rotor voltage");
xlabel("time [s]")
ylabel("u [V]")

# calculate magnetic flux from induction voltage
# Ui=dPhi/dt
# Ui=d(B*A)/dt
# B=1/A*integral(Ui*dt)

offset=mean(vRotor); % the oszi channel has some offset
dt=t(2)-t(1);
# integrator: square integration
numberOfWindings=250;
# iron sheet
dSheet=1;    %Blechdicke
wSheet=14.8; % Blechbreite
d_m=dSheet/1000;
w_m=wSheet/1000;
A=d_m*w_m*numberOfWindings;
B=cumsum(vRotor-offset)/A*dt; 
centerOffset=mean(B)

figure

%plot(B)

BB=B-centerOffset;
plot(iField_A,BB)
hold on
plot(iFieldR_A,BB)
legend("sensor","Rshunt");
grid on
xlabel("current [A]")
ylabel("B [T]")

hold off