% OilDropSpecV4.m
% This program calculates various things from one absorption spectrum.
% Tim Sesterhenn, iposita@gmail.com
% version 1: April 2, 2012
% version 2: April 10, 2012
% version 3: May 25, 2012
% version 4: August 12, 2012

clear a* b* c* d* e* f* g* h* i* j* k* l* m* n* o* p* q* r* s* t* u* v* w* x* y*
% start with everything cleared except for variables starting with z
% (important variables that need to be maintained)
go_out = 0;  % do not write output now
dat = '.dat';

disp(' ')
disp('OilDropSpec version 4, oil drop spectrum analysis program.')
disp(' ')

disp('Give the file path for your input files (enclose in single quotes):')
zdir = input('>> ');
directory = zdir;

% HOW MANY DISCRETE WORDS ARE IN THE SPECIES NAME? (row 12 of the data file)
% e.g.,
% Carduelis tristis AMGO = 3
% Carduelis tristis = 2
% AMGO = 1
sp_words = input...
    ('How many discrete words are in the species description in row 12 of data file? >> ');

while go_out == 0

runtest = exist('z');
if runtest == 0      % if z doesn't exist (the first time run through)
    z = 1;           % make z exist
elseif runtest == 1  % if z does exist (after the first run)
    z = z + 1;       % add 1 to it
end

%-----import data-----

slash = '\';
filename  = ...
    input('Input filename (inside single quotes, no file extension):  >> ');
infile    = [directory slash filename dat];  % makes the whole file path
fid = fopen(infile);  % import data
inputdata = textscan(fid,'%f %c %f %c %f %c %f','headerlines',34);
fclose(fid);
fid = fopen(infile);  % import wavelength range
inputinfo1 = textscan(fid,'%f','headerlines',3);
fclose(fid);
fid = fopen(infile);  % import bird characteristics
inputinfo2 = textscan(fid,'%s','headerlines',12);
fclose(fid);

if sp_words == 1
    fid = fopen(infile);
    inputinfo3 = textscan(fid,'%s',1,'headerlines',11);
    fclose(fid);
    sp1 = inputinfo3{1};
    species = strcat(sp1);
    species = species{1,1};
end
if sp_words == 2
    fid = fopen(infile);
    inputinfo3 = textscan(fid,'%s %s',1,'headerlines',11);
    fclose(fid);
    sp1 = inputinfo3{1};
    sp2 = inputinfo3{2};
    species = strcat(sp1,{' '},sp2);
    species = species{1,1};
end
if sp_words == 3
    fid = fopen(infile);
    inputinfo3 = textscan(fid,'%s %s %s',1,'headerlines',11);
    fclose(fid);
    sp1 = inputinfo3{1};
    sp2 = inputinfo3{2};
    sp3 = inputinfo3{3};
    species = strcat(sp1,{' '},sp2,{' '},sp3);
    species = species{1,1};
end
if sp_words == 4
    fid = fopen(infile);
    inputinfo3 = textscan(fid,'%s %s %s %s',1,'headerlines',11);
    fclose(fid);
    sp1 = inputinfo3{1};
    sp2 = inputinfo3{2};
    sp3 = inputinfo3{3};
    sp4 = inputinfo3{4};
    species = strcat(sp1,{' '},sp2,{' '},sp3,{' '},sp4);
    species = species{1,1};
end
if sp_words == 5
    fid = fopen(infile);
    inputinfo3 = textscan(fid,'%s %s %s %s %s',1,'headerlines',11);
    fclose(fid);
    sp1 = inputinfo3{1};
    sp2 = inputinfo3{2};
    sp3 = inputinfo3{3};
    sp4 = inputinfo3{4};
    sp5 = inputinfo3{5};
    species = strcat(sp1,{' '},sp2,{' '},sp3,{' '},sp4,{' '},sp5);
    species = species{1,1};
end


%-----pick some things out of the data file-----

waves       = inputinfo1{1,1};   % the wavelength range
subjectinfo = inputinfo2{1,1};   % all bird characteristics
nm_min      = waves(1,1);        % minimum wavelength
nm_max      = waves(2,1);        % maximum wavelength
sex         = subjectinfo{1,1};  % sex of subject
birdID      = subjectinfo{2,1};  % band ID
prepnum     = subjectinfo{3,1};  % prep number
subjage     = subjectinfo{4,1};  % subject age
eyeLR       = subjectinfo{5,1};  % which eye (L=left, R=right)
component   = subjectinfo{7,1};  % component
ODtype      = subjectinfo{8,1};  % oil droplet type
curve1      = inputdata{1,5};    % the absorbance data points
nm_range    = [nm_min:1:nm_max]; % makes wavelength array


%-----plot the curve-----

bigcurve = figure('Units','normalized','position',[0.01 0.4 0.5 0.5]);
  % make the plot in the upper-left corner of the monitor
hold on  % so you can plot lines on it later
ticks_here = [nm_min:20:nm_max];
set(gca,'XTick',ticks_here)  % what the x-ticks are
gridx = [nm_min:10:nm_max];  % start making plot gridlines
gridx2= [nm_min+5:10:nm_max];
gridbot = min(curve1)+(0.001*max(curve1));
gridtop = max(curve1)+(0.1*max(curve1));
onearray = ones(1,size(gridx,2));
onearray2= ones(1,size(gridx2,2));
gridybot = gridbot.*onearray;
gridytop = gridtop.*onearray;
gridybot2= gridbot.*onearray2;
gridytop2= gridtop.*onearray2;
for p = 1:1:size(gridx,2)
    plot([gridx(p) gridx(p)],[gridybot(p) gridytop(p)],'color',[0.5 0.5 0.5])
end
for p = 1:1:size(gridx2,2)
    plot([gridx2(p) gridx2(p)],[gridybot2(p) gridytop2(p)],'color',[0.7 0.7 0.7])
end
ysmrng  = [gridbot:0.005:gridtop];
ysmrng2 = [gridbot+0.0025:0.005:gridtop];
for q = 1:1:size(ysmrng,2)
    plot([nm_range(1) nm_range(size(nm_range,2))],...
        [ysmrng(q) ysmrng(q)],'color',[0.5 0.5 0.5])
end
for q = 1:1:size(ysmrng2,2)
    plot([nm_range(1) nm_range(size(nm_range,2))],...
        [ysmrng2(q) ysmrng2(q)],'color',[0.7 0.7 0.7])
end  % end of making gridlines
for t = 1:1:size(ticks_here,2)
    plot([ticks_here(t) ticks_here(t)],[0 0.05*max(curve1)],'-k',...
        'LineWidth',2)
end
plot(nm_range,curve1,'-k')  % plot the basic values


%-----find baseline-----

startbase = ...
    input('Select a wavelength a little past the peak (integer):  >> ');
base_here  = startbase-1;  % subtract 1 to make the loop happy
test_slope = -2;  % initialize it with a loop-friendly value
test_curve = transpose(curve1);  % make absorbance values the right shape
while test_slope < 0.00  % look for where the slope becomes positive
    base_here  = base_here+1;  % see, the loop fixed it
    if base_here == 700  % slope didn't get positive = bad
        disp('Something is wrong with the slope routine. Say that it looks');
        disp('bad and input a new baseline value. Sorry about that.')
        break
    end
    test_range = base_here:1:base_here+15;  % look at start + 15 nm
    test_index = test_range - nm_min;  % convert nm to array indices
    test_eqn   = polyfit(test_index,test_curve(test_index),1);  % line fit
    test_slope = test_eqn(1);  % get the slope of your range of points
end
baseline = curve1(base_here-nm_min);  % convert nm to index
plot([nm_range(1) nm_range(size(nm_range,2))],[baseline baseline],'--r')
is_ok = input('Does the baseline look good? (enter 1 for yes, 0 for no) >> ');
if is_ok == 0
    be_here = input('What absorbance value would be better? >> ');
    plot([nm_range(1) nm_range(size(nm_range,2))],[be_here be_here],'--m')
    is_ok2 = input('Does the magenta line look fine? (1 for yes, 0 for no) >> ');
    if is_ok2 == 0
        be_here = input...
            ('Last chance.  What absorbance value do you REALLY want? >> ');
        plot([nm_range(1) nm_range(size(nm_range,2))],[be_here be_here],'--g')
        disp('The green line is your baseline.  I hope you like it.');
        disp('If you don''t like it, finish this file then redo it.');
    end
    baseline = be_here;  % replace computed baseline with inputted value
end
curve2 = curve1 - baseline;  % curve scaled to baseline


%-----find the peak-----

startpeak = input('At what wavelength should the peak-finding range begin? >> ');
plot([startpeak startpeak],[gridybot(1) gridytop(1)],'color',[0.2 0.56 1])
endpeak   = input('At what wavelength should the peak-finding range end?  >> ');
plot([endpeak endpeak],[gridybot(1) gridytop(1)],'color',[0.2 0.56 1])
peak_range = [startpeak:1:endpeak];  % nm
pk_rng_ind = peak_range - nm_min;    % indices
peak_vals  = test_curve(pk_rng_ind); % absorbances in that range
P_eqn      = polyfit(peak_range,peak_vals,2);  % fit a quadratic to range
peak_curve = (P_eqn(1).*(peak_range.^2))+(P_eqn(2).*peak_range)+P_eqn(3);
plot(peak_range,peak_curve,'color',[0.2 0.56 1])  % quadratic curve
peak = max(peak_curve);  % the top of the quadratic curve

plot([nm_range(1) nm_range(size(nm_range,2))],[peak peak],'--b')
pk_ok = input('Does the peak look good? (enter 1 for yes or 0 for no) >> ');
if pk_ok == 0
    pk_here = input('What absorbance value would be better? >> ');
    plot([nm_range(1) nm_range(size(nm_range,2))],[pk_here pk_here],'--m')
    pk_ok2 = input('Does the magenta line look fine? (1 for yes, 0 for no) >> ');
    if pk_ok2 == 0
        pk_here = input...
            ('One more try.  What is the REAL peak absorbance? >> ');
        plot([nm_range(1) nm_range(size(nm_range,2))],[pk_here pk_here],'--g')
        disp('The green line is your final peak value.  I hope it pleases you.');
        disp('If it displeases you, finish this file then redo it.');
        disp('Better luck next time.');
    end
    peak = pk_here;  % replace computed peak with inputted peak
end
hold off  % stop plotting stuff on the original graph
disp(' ')
disp('Please wait a few seconds while the standardized graph is made.')
pause(5);

%-----plot smaller curve-----

norm_curve = (test_curve-baseline)./(peak-baseline);
littlecurve = figure('Units','normalized','position',[0.01 0.4 0.5 0.5]);
  % make the plot in the upper-left corner of the monitor
hold on  % so you can plot lines on it later
go = startpeak-50;
if go < nm_min+1
    go = nm_min+1;
end
rnd5 = floor(go/5);
go = rnd5 * 5;
nogo = ceil( (base_here+50) /5 );
nogo = nogo * 5;
sm_range = [go:1:nogo];
sm_index = sm_range - nm_min;
if sm_index(1) == 0
    sm_index = sm_index + 1;
end
plot(sm_range,norm_curve(sm_index),'-k','LineWidth',2)  % plot the basic values
set(gca,'XTick',nm_min:20:nm_max)  % where the x-ticks are
set(gca,'XLim',[go-5 nogo+5])
set(gca,'YLim',[-0.1 1.1])
gridx = [go:5:nogo];
gridx2= [go+5:10:nogo];
gridbot = -0.1;
gridtop = 1.1;
onearray = ones(1,size(gridx,2));
onearray2= ones(1,size(gridx2,2));
gridybot = gridbot.*onearray;
gridytop = gridtop.*onearray;
gridybot2 = gridbot.*onearray2;
gridytop2 = gridtop.*onearray2;
for p = 1:1:size(gridx,2)
    plot([gridx(p) gridx(p)],[gridybot(p) gridytop(p)],'color',[0.5 0.5 0.5])
end
for p = 1:1:size(gridx2,2)
    plot([gridx2(p) gridx2(p)],[gridybot2(p) gridytop2(p)],'color',[0.7 0.7 0.7])
end
ysmrng  = [0:0.1:1];
ysmrng2 = [-0.05:0.1:1.05];
for q = 1:1:size(ysmrng,2)
    plot([(go-5) (nogo+5)],[ysmrng(q) ysmrng(q)],'color',[0.5 0.5 0.5])
end
for q = 1:1:size(ysmrng2,2)
    plot([(go-5) (nogo+5)],[ysmrng2(q) ysmrng2(q)],'color',[0.7 0.7 0.7])
end
plot(sm_range,norm_curve(sm_index),'-k','LineWidth',2)
plot([go nogo],[0.0 0.0],'--r','LineWidth',2)
plot([go nogo],[1.0 1.0],'--b','LineWidth',2)
plot([go nogo],[0.5 0.5],'--m','LineWidth',2)


%-----find lambda-mid-----

mid_start = input('Lambda mid is near what wavelength (integer)?  >> ');
Lm_range  = mid_start-10:1:mid_start+10;  % 10 nm to either side of mid
Lmid_X    = Lm_range - nm_min;  % convert range to indices
Lmid_Y    = norm_curve(Lmid_X);  % get the absorbancies for the Lmid range
%fit line
Lmid_eqn   = polyfit(Lm_range,Lmid_Y,1);  % first order (1=linear)
precise_nm = mid_start-10:0.1:mid_start+10;  % nm axis, but more precise
%make line from equation and precise_nm
Lmid_line = (Lmid_eqn(1) .* precise_nm) + Lmid_eqn(2);
plot([precise_nm(1) precise_nm(end)],[Lmid_line(1) Lmid_line(end)],'-c',...
    'LineWidth',3)
minus_mid  = abs(Lmid_line - 0.5);
     % closeness of absorbancies to 0.5
[small,here] = min(minus_mid);
Lmid = precise_nm(here);


%-----calculate stuff-----

Lcut   = (1 - (Lmid_eqn(2))) / Lmid_eqn(1);  % Lcut = (1-intercept)/slope
ntrcpt = Lmid_eqn(2);
Lslope = Lmid_eqn(1);
Bmid   = 0.5/(Lmid-Lcut);
b      = 2.89 * Bmid;
Lzero  = Lcut + (0.37/Bmid);
ypred  = polyval(Lmid_eqn,(Lmid_X + nm_min));  % predicted y's for x-range
SSY    = sum( (Lmid_Y - mean(Lmid_Y) ).^2 );
RSS    = sum((Lmid_Y - ypred).^2);           % residual sum of squares
SSR    = SSY - RSS;
Rsq    = SSR/SSY; % R-squared
disp('R-squared for your line is:');
disp(Rsq);
if Rsq < 0.6
    disp('                                         ');
    disp('****  DANGER!!   DANGER!!  DANGER!!  ****');
    disp('*** THE R-SQUARED VALUE IS BELOW 0.6! ***');
    disp('Rsq = ');
    disp(Rsq);
    disp('****  DANGER!!   DANGER!!  DANGER!!  ****');
    disp('                                         ');
end


%-----make output arrays-----

outnumbers = [Lmid Lcut Lzero b Bmid Rsq ntrcpt Lslope baseline peak];   

if z == 1  % first time through: create the arrays
    zxlsoutspeci = {species};
    zxlsoutsex   = {sex};
    zxlsoutID    = {birdID};
    zxlsoutprep  = {prepnum};
    zxlsoutage   = {subjage};
    zxlsoutLR    = {eyeLR};
    zxlsoutcompo = {component};
    zxlsouttype  = {ODtype};
    zoutcalcs    = outnumbers';
    zfilenames   = {filename};
    zoutnumbers  = norm_curve';
end
if z > 1  % additional times: append new data to the old
    zxlsoutspeci = [zxlsoutspeci {species}];
    zxlsoutsex   = [zxlsoutsex {sex}];
    zxlsoutID    = [zxlsoutID {birdID}];
    zxlsoutprep  = [zxlsoutprep {prepnum}];
    zxlsoutage   = [zxlsoutage {subjage}];
    zxlsoutLR    = [zxlsoutLR {eyeLR}];
    zxlsoutcompo = [zxlsoutcompo {component}];
    zxlsouttype  = [zxlsouttype {ODtype}];
    zoutcalcs    = [zoutcalcs outnumbers'];
    zfilenames   = [zfilenames {filename}];
    zoutnumbers  = [zoutnumbers norm_curve'];
end
    
    
%-----write final Excel output-----

if z == 250
    disp('250 files!?  Wow.  Data will export now so it fits in Excel.');
    go_out = 1;
end
if z < 250
    go_out = input('Are you finished for this session? 1 for yes, 0 for no >> ');
end

close(bigcurve);
close(littlecurve);

end %while loop

if go_out == 1
    disp(' ');
    disp('Writing output file to the same directory as your input files.');
    disp('Please wait until "End of Program" is displayed.');
    disp(' ');
    rightnow = clock;  % set up file name by time of day
    hr = rightnow(4);
    mint = rightnow(5);
    nm_col = nm_range';  % makes nm_range into a column
    goherexls = [directory slash 'ODS_' date '_' num2str(hr) num2str(mint) '.xls'];  % output file
    wordscol = {'Species            ';'Sex                ';...
                'Bird ID            ';'Prep #             ';...
                'Age                ';'Eye                ';...
                'Cellular Component ';'Type               ';...
                'Lambda Mid         ';'Lambda Cut         ';...
                'Lambda 0           ';'b                  ';...
                'Bmid               ';'R Squared          ';...
                'Intercept          ';'Slope              ';...
                'Baseline Absorbance';'Peak Absorbance    ';...
                'Filename           '};
    xlswrite(goherexls,wordscol,'MSPoutput','A1');
    xlswrite(goherexls,nm_col,'MSPoutput','A20');
    xlswrite(goherexls,zxlsoutspeci,'MSPoutput','B1');
    xlswrite(goherexls,zxlsoutsex,'MSPoutput','B2');
    xlswrite(goherexls,zxlsoutID,'MSPoutput','B3');
    xlswrite(goherexls,zxlsoutprep,'MSPoutput','B4');
    xlswrite(goherexls,zxlsoutage,'MSPoutput','B5');
    xlswrite(goherexls,zxlsoutLR,'MSPoutput','B6');
    xlswrite(goherexls,zxlsoutcompo,'MSPoutput','B7');
    xlswrite(goherexls,zxlsouttype,'MSPoutput','B8');
    xlswrite(goherexls,zoutcalcs,'MSPoutput','B9');
    xlswrite(goherexls,zfilenames,'MSPoutput','B19');
    xlswrite(goherexls,zoutnumbers,'MSPoutput','B20');
    if hr < 4
        disp(' ');
        disp('You sure are up late.  Such a hard worker!');
        disp('You should probably go to bed now.  Goodnight.');
        disp(' ');
        disp('End of Program.')
        disp(' ')
    else
        disp(' ');
        disp('Excellent work!  Enjoy the rest of your day!');
        disp(' ');
        disp('End of Program.')
        disp(' ')
    end
    clear z*
end