[Filename,Pathname] = uigetfile('*.txt','Insert Spatial-PL .txt file into program');
cd(Pathname);
SHGdata = importdata(Filename);

%insert the change in HWP angle between scans
ChanDeg = 1;

%get size of raw data array
datasize = size(SHGdata);
dataheight = datasize(1);
datawidth = datasize(2)-1;

%take the first column of the spectrum array
wavelength = SHGdata(1:dataheight,1);

%input wavelength start and end for your signal
WavelengthStart = 400; % you set this
WavelengthStartIndex = 0; %the index associated with WavelengthStart 
WavelengthEnd = 410; % you set this
WavelengthEndIndex = 0; %the index associated with WavelengthEnd 


%find the array index associated with WavelengthStart and WavelengthEnd
z = 0;
while z<=dataheight
    z=z+1;
    if  wavelength(z) >= WavelengthStart
        WavelengthStartIndex = z;
        break;   
    end 
end
while z<=dataheight
    z=z+1;
    if  wavelength(z) >= WavelengthEnd
        WavelengthEndIndex = z;
        break;
    end
end

% uncomment below if you want to adjust for a linear slope to the data
i=1;
adjustment = zeros(dataheight, datawidth+1);
y1=SHGdata(1, 2:datawidth+1);
y2=SHGdata(dataheight-1, 2:datawidth+1);
while i <= dataheight
    adjustment(i, 2:datawidth+1)= ((y2-y1)/(dataheight-1))*(i-1)+y1;
    i=i+1;
end
SHGdata = SHGdata-adjustment;

%uncomment below if you want to subtract the minimum of your data
M = min(SHGdata(WavelengthStartIndex:WavelengthEndIndex,:),[],1);
adjustment = ones(dataheight, datawidth+1);
adjustment(:,1) = 0*adjustment(:,1);
i=1;
while i < datawidth+2
        adjustment(:,i) = adjustment(:,i)*M(i);
        i = i + 1;
end
SHGdata = SHGdata-adjustment;

%take slice of data that includes peak
CutSHGData = SHGdata(WavelengthStartIndex:WavelengthEndIndex, :);

%creates 2xdatawidth array of zeros for output
integratedSHGdata = zeros(2,datawidth); 

%integrate (sum) data and output to array
x=1;
while x <= datawidth
    integratedSHGdata(2,x) = sum(CutSHGData(:,x+1));
    integratedSHGdata(1,x) = 2 * ChanDeg * x; % convert HWP angle to degrees
    x = x + 1;
end

%attempt to fit sine wave
angle = integratedSHGdata(1,:);
y=integratedSHGdata(2,:);

yu = max(y);
yl = min(y);
yr = (yu-yl);             % Range of ‘y’
yz = y-yu+(yr/2);
per = 60   ;              % Estimate period
ym = mean(y);             % Estimate offset

fit = @(b,x)  b(1).*(sin((x - b(2))*2*pi/(per))) + b(3);    % Function to fit
fcn = @(b) sum((fit(b,angle) - y).^2);                      % Least-Squares cost function
s = fminsearch(fcn, [yr;  0;  ym])     ;  % Minimise Least-Squares

figure;
plot(angle, y, angle, fit(s,angle));
%hold on;
%plot(angle, y)
xlabel('angle (degrees')
ylabel('integrated intensity a.u.')
legend('data', 'fit')
s(2);

% that's it!