% This Matlab script quantifies polarotactic responses of tethered flies in
% resonse to a rapidly switching e-vector (90deg back and forth). All .TXT 
% files within the Matlab folder will be processed. The .TXT files are the output
% files of the tracking process using FIJI (macro with instructions can be downloaded at 
% www.flygen.org/skylight-navigation). Please make sure to define 
% the correct data column (Angle1) to be extracted from the .TXT files.
% For each processed file this script will automatically save a .CSV with 
% processed parameters as well as an image for each file containing plots 
% visializing these polarotactic parameters for better overview.
% This script uses the AngleAverage script by Robert Goldstein which can be
% downloaded at:
% https://www.mathworks.com/matlabcentral/fileexchange/26648-angleaverage?focused=5145895&tab=function.



sad = transpose( dir( '*.txt' ) );

    for file = sad
        disp( file.name )
        filename = (file.name)
%% Initialize variables.
delimiter = '\t';

%% Read columns of data as strings:
% For more information, see the TEXTSCAN documentation.
formatSpec = '%*s%s%*s%*s%*s%s%[^\n\r]';

%% Open the text file.
fileID = fopen(filename,'r');

%% Read columns of data according to format string.
% This call is based on the structure of the file used to generate this
% code. If an error occurs for a different file, try regenerating the code
% from the Import Tool.
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter,  'ReturnOnError', false);

%% Close the text file.
fclose(fileID);

%% Convert the contents of columns containing numeric strings to numbers.
% Replace non-numeric strings with NaN.
raw = repmat({''},length(dataArray{1}),length(dataArray)-1);
for col=1:length(dataArray)-1
    raw(1:length(dataArray{col}),col) = dataArray{col};
end
numericData = NaN(size(dataArray{1},1),size(dataArray,2));

% Converts strings in the input cell array to numbers. Replaced non-numeric
% strings with NaN.
rawData = dataArray{2};
for row=1:size(rawData, 1);
    % Create a regular expression to detect and remove non-numeric prefixes and
    % suffixes.
    regexstr = '(?<prefix>.*?)(?<numbers>([-]*(\d+[\,]*)+[\.]{0,1}\d*[eEdD]{0,1}[-+]*\d*[i]{0,1})|([-]*(\d+[\,]*)*[\.]{1,1}\d+[eEdD]{0,1}[-+]*\d*[i]{0,1}))(?<suffix>.*)';
    try
        result = regexp(rawData{row}, regexstr, 'names');
        numbers = result.numbers;
        
        % Detected commas in non-thousand locations.
        invalidThousandsSeparator = false;
        if any(numbers==',');
            thousandsRegExp = '^\d+?(\,\d{3})*\.{0,1}\d*$';
            if isempty(regexp(numbers, thousandsRegExp, 'once'));
                numbers = NaN;
                invalidThousandsSeparator = true;
            end
        end
        % Convert numeric strings to numbers.
        if ~invalidThousandsSeparator;
            numbers = textscan(strrep(numbers, ',', ''), '%f');
            numericData(row, 2) = numbers{1};
            raw{row, 2} = numbers{1};
        end
    catch me
    end
end


%% Split data into numeric and cell columns.
rawNumericColumns = raw(:, 2);
rawCellColumns = raw(:, 1);


%% Replace non-numeric cells with NaN
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),rawNumericColumns); % Find non-numeric cells
rawNumericColumns(R) = {NaN}; % Replace non-numeric cells

%% Allocate imported array to column variable names
Label1 = rawCellColumns(:, 1);
Angle1 = cell2mat(rawNumericColumns(:, 1));


%% Clear temporary variables
clearvars filename delimiter formatSpec fileID dataArray ans raw col numericData rawData row regexstr result numbers invalidThousandsSeparator thousandsRegExp me rawNumericColumns rawCellColumns R;



% Define range of angular data to be analyzed
% -45deg:
a=Angle1([600:1800 4200:5400 7800:9000 11400:12600 15000:16200])
% 45deg:
b=Angle1([2400:3600 6000:7200 9600:10800 13200:14400 16800:17990])

% Average fly heading -45deg
arad=deg2rad(a)
adoublerad=arad*2
adoubleradpiminpi=shiftAnglesFromMinusPIToPI(adoublerad)
meanadoubleradpiminpi = AverageAnglesRadians(adoubleradpiminpi);
meanadoubleradpiminpideg = rad2deg(meanadoubleradpiminpi);
meanadoubleraddeg0360=meanadoubleradpiminpideg + (meanadoubleradpiminpideg < 0)*360
meanadoubleraddeg0180=meanadoubleraddeg0360/2

% Average fly heading 45deg
brad=deg2rad(b)
bdoublerad=brad*2
bdoubleradpiminpi=shiftAnglesFromMinusPIToPI(bdoublerad)
meanbdoubleradpiminpi = AverageAnglesRadians(bdoubleradpiminpi);
meanbdoubleradpiminpideg = rad2deg(meanbdoubleradpiminpi);
meanbdoubleraddeg0360=meanbdoubleradpiminpideg + (meanbdoubleradpiminpideg < 0)*360
meanbdoubleraddeg0180=meanbdoubleraddeg0360/2

% Absolute difference between means of a/b
absdiffbetweenmeans=abs(meanadoubleraddeg0180-meanbdoubleraddeg0180)
absdiffbetweenmeansfinal=min(180-absdiffbetweenmeans, absdiffbetweenmeans)

% Mean difference to e-vector a/b
absdifftoevectora=abs(meanadoubleraddeg0180-45)
absdifftoevectorafinal=min(180-absdifftoevectora, absdifftoevectora)
absdifftoevectorb=abs(meanbdoubleraddeg0180-135)
absdifftoevectorbfinal=min(180-absdifftoevectorb, absdifftoevectorb)
abspreferreddirectiontoevector=(absdifftoevectorafinal+absdifftoevectorbfinal)/2

% Measure how much preferred orientation differs between a and b
absdiffbetweenpreferreddirectiontoevector=abs(absdifftoevectorafinal-absdifftoevectorbfinal)

% Calculate speed data with moving mean
speedarad=deg2rad(Angle1(1:17990))
speedaraddouble=speedarad*2
speedaradunwrap=unwrap(speedaraddouble)
speedaradunwrapdeg=rad2deg(speedaradunwrap)
speedafinal2=diff(speedaradunwrapdeg)
speedafinal3=speedafinal2/2
speedafinal3a=speedafinal3(1:17400)
speedafinaljumpmean=reshape(speedafinal3a,600,[])
yy=sum(speedafinaljumpmean,1)./size(speedafinaljumpmean,1)
speedafinaljumpmeanfinalrow=reshape(repmat(yy, size(speedafinaljumpmean,1),1),1,[])
speedafinaljumpmeanfinal=transpose(speedafinaljumpmeanfinalrow)
meanlastbit=mean(speedafinal3(17401:17989))
meanlastbitlinrow=linspace(meanlastbit,meanlastbit,590)
meanlastbitlin=transpose(meanlastbitlinrow)
speedafinaljumpmeanfinalfinal=[speedafinaljumpmeanfinal; meanlastbitlin]
speedafinaljumpmeanfinalfinala=speedafinaljumpmeanfinalfinal([600:1800 4200:5400 7800:9000 11400:12600 15000:16200])
speedafinaljumpmeanfinalfinalb=speedafinaljumpmeanfinalfinal([2400:3600 6000:7200 9600:10800 13200:14400 16800:17990])
speedAngle1MOVMEAN=movmean(speedafinal3, 10)
speedAngle1absMOVMEAN=abs(speedAngle1MOVMEAN)
speedAngle1absfinalMOVMEAN=min(180-speedAngle1absMOVMEAN, speedAngle1absMOVMEAN)
speedAngle1absfinalsecMOVMEAN=speedAngle1absfinalMOVMEAN*60
speed10befMOVMEAN=speedAngle1absfinalsecMOVMEAN([1200:1800 3000:3600 4800:5400 6600:7200 8400:9000 10200:10800 12000:12600 13800:14400 15600:16200])
speed10aftMOVMEAN=speedAngle1absfinalsecMOVMEAN([1800:2400 3600:4200 5400:6000 7200:7800 9000:9600 10800:11400 12600:13200 14400:15000 16200:16800])
meanspeed10befMOVMEAN=mean(speed10befMOVMEAN)
meanspeed10aftMOVMEAN=mean(speed10aftMOVMEAN)
speed10befchunkMOVMEAN=mean(speedAngle1absfinalsecMOVMEAN([1200:1800; 3000:3600; 4800:5400; 6600:7200; 8400:9000; 10200:10800; 12000:12600; 13800:14400; 15600:16200]))
speed10aftchunkMOVMEAN=mean(speedAngle1absfinalsecMOVMEAN([1800:2400; 3600:4200; 5400:6000; 7200:7800; 9000:9600; 10800:11400; 12600:13200; 14400:15000; 16200:16800]))
speed10chunkfinalMOVMEAN=([speed10befchunkMOVMEAN speed10aftchunkMOVMEAN])
speeddifffinalsecMOVMEAN=meanspeed10aftMOVMEAN-meanspeed10befMOVMEAN


% Number of frames polarotaxis
% Quantifies polarotxis when the angular speed of the fly 
% is smaller than 0.05 deg/frame for 10s
numberframespolarotaxisa=sum(speedafinaljumpmeanfinalfinala>-0.05 & speedafinaljumpmeanfinalfinala<0.05)
number10secondsofpolarotaxisa=numberframespolarotaxisa/600
numberframespolarotaxisb=sum(speedafinaljumpmeanfinalfinalb>-0.05 & speedafinaljumpmeanfinalfinalb<0.05)
number10secondsofpolarotaxisb=numberframespolarotaxisb/600
number10secondsofpolarotaxiswhole=number10secondsofpolarotaxisa+number10secondsofpolarotaxisb
polarotaxisonly=speedafinaljumpmeanfinalfinal>-0.05 & speedafinaljumpmeanfinalfinal<0.05
polarotaxisonly180=polarotaxisonly*180


% Polarotaxis only Section
polarotaxisonlya=polarotaxisonly([600:1800 4200:5400 7800:9000 11400:12600 15000:16200])
polarotaxisonlyb=polarotaxisonly([2400:3600 6000:7200 9600:10800 13200:14400 16800:17990])

Polarotaxisonlya=a(polarotaxisonlya)
aradPolarotaxisonly=deg2rad(Polarotaxisonlya)
adoubleradPolarotaxisonly=aradPolarotaxisonly*2
adoubleradpiminpiPolarotaxisonly=shiftAnglesFromMinusPIToPI(adoubleradPolarotaxisonly)
meanadoubleradpiminpiPolarotaxisonly = AverageAnglesRadians(adoubleradpiminpiPolarotaxisonly);
meanadoubleradpiminpidegPolarotaxisonly = rad2deg(meanadoubleradpiminpiPolarotaxisonly);
meanadoubleraddeg0360Polarotaxisonly=meanadoubleradpiminpidegPolarotaxisonly + (meanadoubleradpiminpidegPolarotaxisonly < 0)*360
meanadoubleraddeg0180Polarotaxisonly=meanadoubleraddeg0360Polarotaxisonly/2

Polarotaxisonlyb=b(polarotaxisonlyb)
bradPolarotaxisonly=deg2rad(Polarotaxisonlyb)
bdoubleradPolarotaxisonly=bradPolarotaxisonly*2
bdoubleradpiminpiPolarotaxisonly=shiftAnglesFromMinusPIToPI(bdoubleradPolarotaxisonly)
meanbdoubleradpiminpiPolarotaxisonly = AverageAnglesRadians(bdoubleradpiminpiPolarotaxisonly);
meanbdoubleradpiminpidegPolarotaxisonly = rad2deg(meanbdoubleradpiminpiPolarotaxisonly);
meanbdoubleraddeg0360Polarotaxisonly=meanbdoubleradpiminpidegPolarotaxisonly + (meanbdoubleradpiminpidegPolarotaxisonly < 0)*360
meanbdoubleraddeg0180Polarotaxisonly=meanbdoubleraddeg0360Polarotaxisonly/2



% Std Deviation from evector within 5min during Poarotaxis for a and b
Polarotaxisonlyarad=deg2rad(Polarotaxisonlya)
Polarotaxisonlyaddouble=Polarotaxisonlyarad*2
Polarotaxisonlyaddoubleradunwrap=unwrap(Polarotaxisonlyaddouble)
Polarotaxisonlyaunwrapdeg=rad2deg(Polarotaxisonlyaddoubleradunwrap)
Polarotaxisonlyaunwrapdeghalf=Polarotaxisonlyaunwrapdeg/2
Polarotaxisonlyaunwrapdeghalfmovmean=movmean(Polarotaxisonlyaunwrapdeghalf,60)
PolarotaxisonlyaStd=std(Polarotaxisonlyaunwrapdeghalf)
PolarotaxisonlyaStd1=std(Polarotaxisonlyaunwrapdeghalfmovmean)

Polarotaxisonlybrad=deg2rad(Polarotaxisonlyb)
Polarotaxisonlybddouble=Polarotaxisonlybrad*2
Polarotaxisonlybddoubleradunwrap=unwrap(Polarotaxisonlybddouble)
Polarotaxisonlybunwrapdeg=rad2deg(Polarotaxisonlybddoubleradunwrap)
Polarotaxisonlybunwrapdeghalf=Polarotaxisonlybunwrapdeg/2
Polarotaxisonlybunwrapdeghalfmovmean=movmean(Polarotaxisonlybunwrapdeghalf,60)
PolarotaxisonlybStd=std(Polarotaxisonlybunwrapdeghalf)
PolarotaxisonlybStd1=std(Polarotaxisonlybunwrapdeghalfmovmean)


% Absolute difference between means of a/b
absdiffbetweenmeansPolarotaxisonly=abs(meanadoubleraddeg0180Polarotaxisonly-meanbdoubleraddeg0180Polarotaxisonly)
absdiffbetweenmeansfinalPolarotaxisonly=min(180-absdiffbetweenmeansPolarotaxisonly, absdiffbetweenmeansPolarotaxisonly)
absdifftoevectoraPolarotaxisonly=abs(meanadoubleraddeg0180Polarotaxisonly-45)
absdifftoevectorafinalPolarotaxisonly=min(180-absdifftoevectoraPolarotaxisonly, absdifftoevectoraPolarotaxisonly)
absdifftoevectorbPolarotaxisonly=abs(meanbdoubleraddeg0180Polarotaxisonly-135)
absdifftoevectorbfinalPolarotaxisonly=min(180-absdifftoevectorbPolarotaxisonly, absdifftoevectorbPolarotaxisonly)
abspreferreddirectiontoevectorPolarotaxisonly=(absdifftoevectorafinalPolarotaxisonly+absdifftoevectorbfinalPolarotaxisonly)/2


% Measure for how much preferred orientation changes between a and b
absdiffbetweenpreferreddirectiontoevectorPolarotaxisonly=abs(absdifftoevectorafinalPolarotaxisonly-absdifftoevectorbfinalPolarotaxisonly)


% Create figures
figure('DefaultAxesFontSize',5)
ax2=subplot(5,3,[1 3])
plot(Angle1)
ax = gca
ax.XTick = [1 1800 3600 5400 7200 9000 10800 12600 14400 16200 18000]
ax.YTick = [0 45 90 135 180]
ax.FontSize = 5
ax.XAxis.Exponent = 0
ax.XAxis.TickLabelFormat = '%6.0f'
axis tight
ylabel('Heading');
ylim([0 180]);
xlim([0 18000])
grid('on');
hold on 
area(polarotaxisonly180,'FaceAlpha',0.1,'EdgeAlpha',0)
hold off


subplot(5,3,4)
Angle1a=Angle1([1:1800 3601:5400 7201:9000 10801:12600 14401:16200])
Angle1b=Angle1([1801:3600 5401:7200 9001:10800 12601:14400 16201:17990])
Angle1rad=deg2rad(Angle1)
Angle1arad=deg2rad(Angle1a)
Angle1brad=deg2rad(Angle1b)
rose(Angle1rad,180)


subplot(5,3,5)
rose(Angle1arad,180)


subplot (5,3,6)
rose(Angle1brad,180)


subplot(5,3,[7])
plot(speed10chunkfinalMOVMEAN)
ax = gca
grid on
ax.FontSize = 3
ax.XAxis.Exponent = 0
ax.XAxis.TickLabelFormat = '%6.0f'
axis([0 1200 0 300])
ax.XTick = [0 600 1200]


subplot(5,3,[8 9])


subplot (5,3,10)
bar([meanadoubleraddeg0180 meanbdoubleraddeg0180 absdiffbetweenmeansfinal])
ax = gca
ax.XTickLabel = ({'a','b','diff'})
ax.FontSize = 5
ax.XAxis.Exponent = 0
ax.XAxis.TickLabelFormat = '%6.0f'
ylabel('Heading Mean')
ylim([0 180]);
ax.YTick = [0 45 90 135 180]
grid('on');


annotation('textbox',...
    [0.408333333333333 0.0888888888888952 0.377380952380952 0.312698412698413],...
    'String',{['Mean Difference= ' num2str(absdiffbetweenmeansfinal) '°'],['Speed Difference 10sec bef/aft= ' num2str(speeddifffinalsecMOVMEAN) '°/sec'],['Diff from Evector a= ' num2str(absdifftoevectorafinal) '°'],['Diff from Evector b= ' num2str(absdifftoevectorbfinal) '°'],['Diff between a/b Diff from Evector (a-b)= ' num2str(absdiffbetweenpreferreddirectiontoevector) '°'],['Mean diff from evector a/b= ' num2str(abspreferreddirectiontoevector) '°'],['Number of 10s Polarotaxis a= ' num2str(number10secondsofpolarotaxisa) '/10'],['Number of 10s Polarotaxis b= ' num2str(number10secondsofpolarotaxisb) '/10'],['Number of 10s Polarotaxis Whole= ' num2str(number10secondsofpolarotaxiswhole) '/20'],['Diff from Evector a Polarotaxisonly= ' num2str(absdifftoevectorafinalPolarotaxisonly) '°'],['Diff from Evector b Polarotaxisonly= ' num2str(absdifftoevectorbfinalPolarotaxisonly) '°'],['Diff between a/b Diff from Evector (a-b) Polarotaxisonly= ' num2str(absdiffbetweenpreferreddirectiontoevectorPolarotaxisonly) '°'],['Mean diff from evector a/b Polarotaxisonly= ' num2str(abspreferreddirectiontoevector) '°'],['Std Deviation a= ' num2str(PolarotaxisonlyaStd1) '°'],['Std Deviation b= ' num2str(PolarotaxisonlybStd1) '°']},...
  'FontSize',5,...
'FitBoxToText','on');


% Create matrix containing values you would like to save
CompleteVariables=[absdiffbetweenmeansfinal;speeddifffinalsecMOVMEAN; absdifftoevectorafinal;absdifftoevectorbfinal; absdiffbetweenpreferreddirectiontoevector; abspreferreddirectiontoevector; number10secondsofpolarotaxisa; number10secondsofpolarotaxisb; number10secondsofpolarotaxiswhole; absdifftoevectorafinalPolarotaxisonly; absdifftoevectorbfinalPolarotaxisonly; absdiffbetweenpreferreddirectiontoevectorPolarotaxisonly; abspreferreddirectiontoevector; PolarotaxisonlyaStd1; PolarotaxisonlybStd1]
CompleteVariablesRow=CompleteVariables.';

% Save plot as .JPG
filenamestring = regexprep(Label1{3,1},'.avi','.jpg');
filename1 = [filenamestring];
filename2 = sprintf('%s.csv',[filenamestring]);
saveas(gcf, filename1, 'jpg');

% Save results as .CSV
csvwrite(filename2,CompleteVariablesRow)
    end