function [ ] = CompareWDMmodes( )
%Compares the frequency response of outcoupling for Red green and blue in
%Unpack updated to give angular deflection 10/2/14
close all;

% User Defined Variables
series='A';
sample=16;
modes=[1];
transducer='T1';

modenum=length(modes);
figposvec=FigPosition();
figure('position',figposvec);
hold on;
for mindex=1:modenum
    modeint=modes(mindex);
    mode=strcat('M',num2str(modeint));
    RData=UnPackWDM(series, sample, 'Red', mode, transducer);
    GData=UnPackWDM(series, sample, 'Green', mode, transducer);
    BData=UnPackWDM(series, sample, 'Blue', mode, transducer);
    
    %frequency graph
    
    windownum=mindex+(mindex-1); %1,3,5,7--the positions of the freq plots in subplot
    %figure;
    %hold on;
    
    Rpower=VectorNorm(RData.fprojsmooth);
    Gpower=VectorNorm(GData.fprojsmooth);
    Bpower=VectorNorm(BData.fprojsmooth);
   
    subplot(modenum,2,windownum);
    title(strcat('Mode ',num2str(mindex)));
    xlabel('Frequency (MHz)');
    hold on;
    plot(RData.fvals,Rpower,'r');
    plot(GData.fvals,Gpower,'g');
    plot(BData.fvals,Bpower,'b');
    
    %position graph
    subplot(modenum,2,windownum+1);%2,4,6,8 positions of position graphs in subplot
    title(strcat('Mode ',num2str(mindex)));
    xlabel('Angle (deg)');
    axis([-45 0 0 1]);%see of this works
    hold on;
    Rpower=VectorNorm(RData.xproj);
    Gpower=VectorNorm(GData.xproj);
    Bpower=VectorNorm(BData.xproj);
    
    
    plot(RData.tvals,Rpower,'r');
    plot(GData.tvals,Gpower,'g');
    plot(BData.tvals,Bpower,'b');
    
    
    
end
samplename=strcat(series,num2str(sample));
suptitle(strcat('Sample ',samplename));

% write plot to file in sample folder
cd(strcat(series,num2str(sample)));
saveas(gcf,strcat('Sample_',samplename,'_',transducer,'_Comparison_Data.jpeg'));
cd ../
end

function [normv] = VectorNorm(v)
%normalize vector
Denom=max(v);
normv=v./Denom;
out=normv;
end


%*********************


function out = UnPackWDM(series,sample,color,mode,transducer)
%WDM comes in three columns power, freq and position
%   this function reads the file and makes a power matrix for surf

device=strcat(series,num2str(sample));
experiment=strcat(device,'_',color,'_',mode,'_',transducer);
currentpath=strcat(device,'\',experiment,'\');

%currentpath='W8\W8_Blue_M1_T3\';
querypath=strcat('\',device,'\',experiment,'\','*data.csv');
% csvstruct=dir(fullfile(cd,'\W8\W8_Blue_M1_T3\*data.csv'));
csvstruct=dir(fullfile(cd,querypath));
csvfile=strcat(currentpath,csvstruct.name);

isthereacsvfile=FileExistance(csvfile);
if isthereacsvfile == 2
    M=csvread(csvfile,1,0);
    %position will be cols, freq will be rows
    coldim=sum(M(:,3)==M(1,3));%number of positions tested by seeing how many position entries are the same as the first
    %this is also equal to the number of frequencies tested
    
    %make matrix of power values
    powervals=M(:,1);
    P=reshape(powervals,coldim,[]);
    P=double(P);
    
    %filter
    isthereafilter=FileExistance(strcat(currentpath,'filter.bmp'));
    if isthereafilter > 0
        filterpath=strcat(currentpath,'filter.bmp');
        P=PlotFilter(P,filterpath);
    end
    
    %eliminate negative numbers
    P(P<0)=0;
    freqvals=M(1:coldim,2);
    posvals=-M(1:coldim:end,3);%negative cause stage moves in negative direction
else
    P=csvread('NullP.csv');
    freqvals=csvread('fvals.csv');
    posvals=csvread('xvals.csv');
    powervals=[];
end


csvconfigstruct=dir(fullfile(cd,strcat('\',currentpath,'*config.csv')));

pcsvfile=strcat(currentpath,csvconfigstruct.name);
isthereapcsvfile=FileExistance(pcsvfile);
if isthereapcsvfile==2
    fileID=fopen(pcsvfile); % Open file
    absolutepos = textscan(fileID,'%s','delimiter',','); % parse file by commas
    pcbpos=str2num(absolutepos{1}{9}(findstr(absolutepos{1}{9},': ')+2:length(absolutepos{1}{9}))); % pull out the offset from 0 (mm)
    guidedpos=str2num(absolutepos{1}{7}(findstr(absolutepos{1}{7},': ')+2:length(absolutepos{1}{7}))); % pull out the distance to detector (mm)
    thetavals=zeros(size(posvals));
    for pindex=1:length(posvals)
        y=30;%distance in mm from sample to sensor
        x=posvals(pindex)-pcbpos;%this is the x distance to the sensor from the pcb
        thetavals(pindex)=atand(x/y);
    end
    fclose(fileID);
else
    disp('cannot find pcsv file')
    thetavals=zeros(size(posvals));
end

% pcsvfile=strcat(currentpath,'position.csv');
% isthereapcsvfile=FileExistance(pcsvfile);
% if isthereapcsvfile==2
%     disp('found pcsv file')
%     absolutepos=csvread(strcat(currentpath,'position.csv'),1,0);
%     pcbpos=absolutepos(1);
%     guidedpos=absolutepos(2);
%     thetavals=zeros(size(posvals));
%     for pindex=1:length(posvals)
%         y=30;%distance in mm from sample to sensor
%         x=posvals(pindex)-pcbpos;%this is the x distance to the sensor from the pcb
%         thetavals(pindex)=atand(x/y);
%     end
%     
% else
%     disp('cannot find pcsv file')
%     thetavals=zeros(size(posvals));
% end


fproj=max(P,[],2);
xproj=max(P,[],1);
fprojsmooth=smooth(fproj,'moving');



% csvwrite('NullP.csv',zeros(size(P)));
% csvwrite('fvals.csv',freqvals);
% csvwrite('xvals.csv',posvals);

%generate plots, if they don't exist
isthereplots=FileExistance(strcat(currentpath,'Master.png'));
istheredata=isthereacsvfile;
if isthereplots < 1 && istheredata == 2
    WriteFilterBlank(P,currentpath);
    PowerSurface(P,posvals,thetavals,freqvals, strcat(currentpath,'Surf'));
    SpaceProjection(P,posvals,thetavals,strcat(currentpath,'Xproj'));
    FreqProjection(P,freqvals,strcat(currentpath,'Fproj'));
    FreqProjectionSmooth(P,freqvals, strcat(currentpath,'Fproj'));
    MasterPlot(P,posvals,freqvals, strcat(currentpath,'Master'));
   % DynamicFig(P,posvals,freqvals, strcat(currentpath,'Master'));
end

out=struct('P',P,...
    'fvals',freqvals,...
    'xvals',posvals,...
    'tvals',thetavals,...
    'pvals',powervals,...
    'fproj',fproj,...
    'fprojsmooth',fprojsmooth,...
    'xproj',xproj,...
    'currentpath',currentpath);


end


function [] = WriteFilterBlank(P,currentpath)
%function wraps exist and makes sure it only checks the current dir.
%normalize to 256
N=max(max(P));
P=(P./N).*256;
isthereafolder=FileExistance(currentpath);
if isthereafolder == 7
    imwrite(P,strcat(currentpath,'P_image.bmp'));
end
end

function ExistCheck = FileExistance( filename )
%function wraps exist and makes sure it only checks the current dir.
ExistCheck=exist(fullfile(cd,filename),'file');

end


function [P] = PlotFilter(P,filterpath)
%filters data

filter=imread(filterpath);%grayscale
filter=filter./255;%grayscale to binary
filtersize=size(filter);
if filter(3)==3% this means the filter has three channels
    disp('filter image is rgb mode, change to grayscale')
end
%change both to double precision
filter=double(filter);
P=P.*filter;%filter out single mode data

end

function [figposvec] = FigPosition()
%returns position vector for figure
scnsize = get(0,'ScreenSize');
scnwidth=scnsize(3);
scnheight=scnsize(4);
BottomLeftCornerX=0;%scnwidth/8;
BottomLeftCornerY=0;%scnwidth/8;
FigWidth=scnwidth/2;
FigHeight=scnheight/.6;
figposvec=[BottomLeftCornerX,BottomLeftCornerY,FigWidth,FigHeight];

end

function [] = DynamicFig(P,posvals,freqvals, filename)
%returns position vector for figure
h=figure('visible','on');
denom=max(max(P));
P=P./denom;%normalize

[x,y] = ndgrid(posvals,freqvals);
F = griddedInterpolant(x, y, P, 'cubic');

newx=linspace(posvals(1),posvals(end),500);
newy=linspace(freqvals(1),freqvals(end),500);

[Xq,Yq] = ndgrid(newx,newy);
newP=F(Xq,Yq);
surf(Xq,Yq,newP,'EdgeColor','none','LineStyle','none','FaceLighting','phong');
%surf(posvals,freqvals,P,'EdgeColor','none','LineStyle','none','FaceLighting','phong');
%colormap('default');
ylabel('Frequency (MHz)')
xlabel('Angle (deg)')
zlabel('Power (W)')
saveas(h,filename);

%close(h);
end



function [] = MasterPlot(P,posvals,freqvals, filename)
%makes subplot surface plot

scnsize = get(0,'ScreenSize');
scnwidth=scnsize(3);
BottomLeftCornerX=scnwidth/8;
BottomLeftCornerY=scnwidth/8;
FigWidth=scnwidth/4;
FigHeight=scnwidth/3;

h=figure('visible','off','Position',[BottomLeftCornerX,BottomLeftCornerY,FigHeight,FigWidth]);

fproj=max(P,[],2);
xproj=max(P,[],1);
%smoothedfreq=smooth(fproj,'moving');

subplot(2,2,1)

surf(posvals,freqvals,P,'EdgeColor','none','LineStyle','none','FaceLighting','phong');
view([90 90]);
%colormap('hot');
ylabel('Frequency (MHz)');
xlabel('Angle(deg)');

subplot(2,2,2)
plot(posvals,xproj);
xlabel('Angle (deg)')
ylabel('Power (W)')

subplot(2,2,3)
plot(freqvals,fproj);
xlabel('Frequency (MHz)')
ylabel('Power (W)')

subplot(2,2,4)
S=surf(posvals,freqvals,P,'EdgeColor','none','LineStyle','none','FaceLighting','phong');
%colormap('default');
ylabel('Frequency (MHz)')
xlabel('Angle (deg)')
zlabel('Power (W)')

saveas(h,filename,'png');



end


function [] = PowerSurface(P,posvals,thetavals,freqvals, filename)
%makes surface plot
h=figure('visible','off');
colormap('default');
surf(posvals,freqvals,P,'EdgeColor','none','LineStyle','none','FaceLighting','phong');
xlabel('Frequency (MHz)')
ylabel('Position (mm)')
zlabel('Power (W)')
saveas(h,filename,'png');

h=figure('visible','off');
colormap('default');
surf(thetavals,freqvals,P,'EdgeColor','none','LineStyle','none','FaceLighting','phong');
xlabel('Frequency (MHz)')
ylabel('Angle (deg)')
zlabel('Power (W)')
saveas(h,strcat(filename,'angle'),'png');
end

function [] = SpaceProjection(P,posvals,thetavals, filename)
%takes space projeciton of power matrix and saves image of plot to file
xproj=max(P,[],1);

h=figure('visible','off');
plot(thetavals,xproj);
xlabel('Angle (deg)')
ylabel('Power (W)')
saveas(h,strcat(filename,'angle'),'png');

h=figure('visible','off');
plot(posvals,xproj);
xlabel('Position (mm)')
ylabel('Power (W)')
saveas(h,strcat(filename,'space'),'png');

end

function [] = FreqProjection(P,freqvals, filename)
%takes freq projeciton (raw) of power matrix and saves image of plot to file
fproj=max(P,[],2);

h=figure('visible','off');
plot(freqvals,fproj);
xlabel('Frequency (MHz)')
ylabel('Power (W)')
saveas(h,filename,'png');

end

function [] = FreqProjectionSmooth(P,freqvals, filename)
%takes freq projeciton (smoothed) of power matrix and saves image of plot to file
fproj=max(P,[],2);
smoothedfilename=strcat(filename, 'smoothed');
smoothedfreq=smooth(fproj,'moving');

h=figure('visible','off');
plot(freqvals,smoothedfreq);
xlabel('Frequency (MHz)')
ylabel('Power (W)')
saveas(h,smoothedfilename,'png');

end


function out = MakeHot(color)
MR=[0,0; 
    0.02,0.3; %this is the important extra point
    0.3,1;
    1,1];
MG=[0,0;
    0.3,0; 
    0.7,1;
    1,1];
MB=[0,0; 
    0.7,0;
    1,1];

if strcmp(color,'r') || strcmp(color,'Red')
hotR = colormapRGBmatrices(500,MR,MG,MB);%redhot
out=hotR;
end

if strcmp(color,'g') || strcmp(color,'Green')
hotG = colormapRGBmatrices(500,MG,MR,MB);%greenhot
out=hotG;
end

if strcmp(color,'b') || strcmp(color,'Blue')
hotB = colormapRGBmatrices(500,MB,MG,MR);%bluehot
out=hotB;
end

end

function mymap = colormapRGBmatrices( N, rm, gm, bm)
%code by cresspahl from http://cresspahl.blogspot.com/2012/03/expanded-control-of-octaves-colormap.html
  x = linspace(0,1, N);
  rv = interp1( rm(:,1), rm(:,2), x);
  gv = interp1( gm(:,1), gm(:,2), x);
  mv = interp1( bm(:,1), bm(:,2), x);
  mymap = [ rv', gv', mv'];
  %exclude invalid values that could appear
  mymap( isnan(mymap) ) = 0;
  mymap( (mymap>1) ) = 1;
  mymap( (mymap<0) ) = 0;
end
