

%% Data aquisition from Altera DE2 board through Serial port
% Tomography 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% This program is written by Bebeh and Kiko on July 2017.                 %
% The goal is to take data from Altera DE2 board through a serial port.   %
% ALTERA sends out 8 numbers that are 32-bit numbers. They correspond to 8%
% counters which are A, B, A' and B' singles as well as AB, A'B, AB' and  %
% A'B' coincidences respectively.                                         %
%_________________________________________________________________________%
%                                                                         %
%*************************************************************************%
%_________________________________________________________________________%
%                                                                         %
% Each 32-bit number is divided into 7-bit partitions which makes it 5    %
% bytes of 7-bit figures per number. Transmitting through the serial port %
% each number becomes 5 8-bit bytes + a termination byte. Eventually      %
% reading the data with MatLab using 19200 baudrate we receive            %
% 8 (# of countes)*5(bytes per number)+1= 41 numbers in the array.        %
% In order to calculate the actual number for a counter, say A, we mu     %
% multiply the corresponding figures by 2^0, 2^7, 2^14, 2^21 and 2^28.    %
% We save the numbers and plot them.                                      %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;
clear;
close all;
format compact


%% Loop
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prompt = {'Enter the COM# of the Counter:','time interval','Name of file = '};
dlg_title = 'Tomography general parameters';
defaultans = {'COM2','1','Tom'};%[1 length(dlg_title)+10],[1 5;1 5;1 5;1 30]
userinput = inputdlg(prompt,dlg_title,[1 length(dlg_title)+30],defaultans);
counterportnum = userinput{1};
scounter = serial(counterportnum,'BaudRate',19200,'DataBits',8,'StopBits',1,'Parity','none');
timeinterval = str2double(userinput{2}); % time interval for each measurement in seconds
nams = userinput{3};
datanam = nams;  %'TM11805Normal3632'; % data name
dateexp = date; % date
numofstates = 16; % # of states to do the tomography 
numofmeasurements = 1; % # of measurements per each state 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
loop = numofmeasurements; 
timeaxislimit = numofstates*numofmeasurements*timeinterval;
deltat = 40e-9; % pulse width to calculate accidental coincidences
count=1; pausetime=0; time=zeros(loop);statepause=0;
clockt = fix(clock); % saving the initial date/time into a matrix
myDatadimension=41*(timeinterval*10)+40; % timeinterval*10=timeinterval in seconds 
cleandatadimension=41*(timeinterval*10);
myData=zeros(myDatadimension,1);
resultsmatrix=zeros(numofmeasurements,8);
% erasingmatrix(1:(numofmeasurements+1)*numofstates+3,1:8)="";% matrix defined to erase the excel sheet
Sheet1=strcat('Data points',num2str(clockt(1,4:6)));
Sheet2=strcat('Average results',num2str(clockt(1,4:6)));
% Start the excel file to write the gradual results
xlrange1='A1';
warning('off','MATLAB:xlswrite:AddSheet');
% to suppress the warning when the sheet name is not in excel file.
Header1={'A','B','AB','Accidenals'}; 
xlswrite('Tomographyresults.xlsx',Header1,Sheet1,xlrange1)
% End of writing the header for the "Gradual Results" sheet in excel file.
%-------------------------------------
% Start the excel file to write the "Total Results" sheet

warning('off','MATLAB:xlswrite:AddSheet');
% to suppress the warning when the sheet name is not in excel file.
Header2={'A','B','AB','Accidentals'}; 
xl2range1='A1';
xlswrite('Tomographyresults.xlsx',Header2,Sheet2,xl2range1)
% End of writing the header for the "Gradual Results" sheet in excel file.


%% Figure adjustments
screensize = get( groot, 'Screensize' ); %getting screen size
position=[1 screensize(1,4)/2-100 screensize(1,3) screensize(1,4)/2];
f1=figure('Name','Tomography Data Recording','numbertitle','off','Position',screensize,'color',[0.7 0.7 0.7]);
% hold on
%_____________________________________
% Axes Properties
% axes('position',[left bottom width height])

% Axis for Header
axheader=axes('position',[0.45 0.7 0.1 0.05],'visible','off');
axheader.Title.Visible = 'on';
set(get(gca,'title'),'color','w','background','b')% figure header text:white, background:blue
%_____________________________________
% Axes for plots
ax1 = axes('position',[0.08 0.2 0.25 0.3]); % Axies 1 position in the figure
set(get(ax1,'title'),'background','y')
ax1.XLim = [0 numofstates*timeinterval+1];
% ax1.XLabel.String  = 'time (sec)';
% ax1.XLabel.String  ='HH HV VV VH RH RV DV DH DR DD RD HD VD VL HL RL DA';
ax1.XLabel.FontWeight = 'bold';
ax1.XLabel.FontSize = 12; ax1.XLabel.FontName = 'TimesNewRoman';
% set(ax1,'XTickL',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17],{'HH','HV','VV','VH','RH','RV','DV','DH','DR','DD','RD','HD','VD','VL','HL','RL','DA'})
set(ax1,'XTick',1:1*timeinterval:numofstates*timeinterval)
set(ax1,'XTicklabel',{'HH','HV','VV','VH','RH','RV','DV','DH','DR','DD','RD','HD','VD','VL','HL','RL','DA'});
ax1.YLabel.String  = 'singles A';
ax1.YLabel.FontWeight = 'bold';
ax1.YLabel.FontSize = 25;ax1.YLabel.FontName = 'TimesNewRoman';
ax1.YLim = [0 inf];
grid(ax1,'on');
hold(ax1,'on')
%____________________________________

ax2 = axes('position',[0.40 0.2 0.25 0.3]); % Axies 2 position in the figure
set(get(ax2,'title'),'background','y')
ax2.XLim=[0 numofstates*timeinterval+1];
% ax2.XLabel.String  = 'time (sec)';
% ax2.XLabel.String  ='HH HV VV VH RH RV DV DH DR DD RD HD VD VL HL RL DA';
ax2.XLabel.FontWeight = 'bold';
ax2.XLabel.FontSize = 12;ax2.XLabel.FontName = 'TimesNewRoman';
set(ax2,'XTick',1:1*timeinterval:numofstates*timeinterval)
set(ax2,'XTicklabel',{'HH','HV','VV','VH','RH','RV','DV','DH','DR','DD','RD','HD','VD','VL','HL','RL','DA'})
ax2.YLabel.String  = 'singles B';ax2.YLabel.FontWeight = 'bold';
ax2.YLabel.FontSize = 25;ax2.YLabel.FontName = 'TimesNewRoman';
ax2.YLim = [0 inf];
grid(ax2,'on')
hold(ax2,'on')
%____________________________________

ax3=axes('position',[0.72 0.2 0.25 0.3]); % Axies 3 position in the figure
set(get(ax3,'title'),'background','y')
ax3.XLim=[0 numofstates*timeinterval+1];
% ax3.XLabel.String  = 'time (sec)';
% ax3.XLabel.String  ='HH HV VV VH RH RV DV DH DR DD RD HD VD VL HL RL DA';
ax3.XLabel.FontWeight = 'bold';
ax3.XLabel.FontSize = 12;ax3.XLabel.FontName = 'TimesNewRoman';
set(ax3,'XTick',1:1*timeinterval:numofstates*timeinterval)
set(ax3,'XTicklabel',{'HH','HV','VV','VH','RH','RV','DV','DH','DR','DD','RD','HD','VD','VL','HL','RL','DA'})
ax3.YLabel.String  = 'Coinc. AB';ax3.YLabel.FontWeight = 'bold';
ax3.YLabel.FontSize = 25;ax3.YLabel.FontName = 'TimesNewRoman';
ax3.YLim=[0 inf];
grid(ax3,'on')
hold(ax3,'on')
%______________________________________________________________

for stateindexi=1:numofstates
prompt = {'Set the wave plates and then click OK'};
dlg_title = 'Wave plate Orientation';
userinput = inputdlg(prompt,dlg_title,[1 length(dlg_title)+30]);

%-----------------------------------------------------------   
% header for each state measurement results in excel file "Gradual Results"
countt=num2str((stateindexi-1)*(numofmeasurements+1)+2);% to go two lines further (count+1)in excel (because of the header)
xlrange2=strcat('A',countt);
stateindexit=num2str(stateindexi);
stateheader={'state #',stateindexit};
xlswrite('Tomographyresults.xlsx',stateheader,Sheet1,xlrange2);
% header for each state measurement results in excel file "Total Results"
countt=num2str((stateindexi-1)*2+2);
xl2range2=strcat('A',countt);
xlswrite('Tomographyresults.xlsx',stateheader,Sheet2,xl2range2);

% %______________________________________________________________
%% The Loop
% defining object s for serial instrument. BaudRate=19200 bps, DataBits=8
% StopBits=1, Parity=none. 
% The COM port is determined by the Device Manager in Windows.
 
fopen(scounter);  % open the serial port before the inner loop begins.
%****************************************************************************

while ~isequal(count,loop+1)
    
    % Serial data accessing 
    
    for i=1:timeinterval
         myData1 = fread(scounter,512,'uint8'); % reading # of bytes
         myData(1,(i-1)*512+1:i*512) = myData1';
    end 
    
    % finding terminationbyte if the 41th element of myData is not 255
    tbi=0;
    if myData(1,41)~=255
        for i=1:40
           if myData(1,i)==255
               terminationbyteindex=i;
           end
        end 
        tbi=terminationbyteindex;
    end 
    % saving myData portion into cleandata so the array starts with A
    % that is right after the first termination byte (255)
    
    cleandata=myData(1,tbi+1:tbi+cleandatadimension);
    numofcounts=zeros(1,8);
    CD=cleandata; % just to use a shorthand notation CD
    kmax=timeinterval*10; % loop repetation numner for each counter
    L=0;j=0;
    for i=1:8
        j=0;
        for k=1:kmax         
            numofcounts(1,i)=numofcounts(1,i)+CD(1,1+j+L)+2^7*CD(1,2+j+L)+2^14*CD(1,3+j+L)+2^21*CD(1,4+j+L)+2^28*CD(1,5+j+L);
            j=j+41; % the corresponding figure after a tenth of a second
        end
        L=L+5; % next counter partition starts at the next 5th byte
    end 
    numofcountsA=numofcounts(1,1);
    n1(stateindexi)=numofcountsA;
    numofcountsB=numofcounts(1,2);
    n2(stateindexi)=numofcountsB;
    numofcountsAprime=numofcounts(1,3);
    numofcountsBprime=numofcounts(1,4);
    numofcountsAB=numofcounts(1,5);
    n(stateindexi)=numofcountsAB;
    numofcountsAprimeB=numofcounts(1,6);
    numofcountsABprime=numofcounts(1,7);
    numofcountsAprimeBprime=numofcounts(1,8);
    
    %% Plotting the data points on different subplots  
    time(count) = (stateindexi-1)*numofmeasurements*timeinterval+count*timeinterval; % x-axis (time) in seconds   
    % plotting A    
      plot(ax1,time(count),numofcounts(1,1),'. b','MarkerSize',20)
    %----------------------------------------------------------------------
    % plotting B
      plot(ax2,time(count),numofcounts(1,2),'. b','MarkerSize',20)
    %----------------------------------------------------------------------
    % plotting AB
      plot(ax3,time(count),numofcounts(1,5),'. b','MarkerSize',20)
    %----------------------------------------------------------------------
    % Drawing Y-data of the three plots (A, B, AB) at the same time
    descriptionA = num2str(numofcounts(1,1));    
    title(ax1,descriptionA,'FontWeight','bold','FontSize',40,'FontName','Times New Roman');
    
    descriptionB = num2str(numofcounts(1,2));
    title(ax2,descriptionB,'FontWeight','bold','FontSize',40,'FontName','Times New Roman');
    
    descriptionAB = num2str(numofcounts(1,5));
    title(ax3,descriptionAB,'FontWeight','bold','FontSize',40,'FontName','Times New Roman');
    
    descr = ['Tomography: ','State # ',num2str(stateindexi),'/',num2str(numofstates),', measurement #',num2str(count),'/',num2str(numofmeasurements)];
    title(axheader,descr,'FontWeight','bold','FontSize',30,'FontName','Times New Roman')
%     text(axheader,0,0,descr,'FontWeight','bold','FontSize',30,'FontName','Times New Roman')
    drawnow

    %----------------------------------------------------------------------
    
    
    %% Storing Data in the resultsmatrix and finally in xlsx file
    
    resultsmatrix(count,1)=numofcounts(1,1);
    resultsmatrix(count,2)=numofcounts(1,2);
    resultsmatrix(count,3)=numofcounts(1,3);
    resultsmatrix(count,4)=numofcounts(1,4);
    resultsmatrix(count,5)=numofcounts(1,5);
    resultsmatrix(count,6)=numofcounts(1,6);
    resultsmatrix(count,7)=numofcounts(1,7);
    resultsmatrix(count,8)=numofcounts(1,8);
%________________________________________________________________________    
    % writing results gradually into the "Gradual Results' sheet in excel
    % file
    warning('off','MATLAB:xlswrite:AddSheet');
    % to suppress the warning when the sheet name is not in excel file.
    countt=num2str((stateindexi-1)*(numofmeasurements+1)+count+2);% to go two lines further (count+1)in excel (because of the header)
    xlrange2=strcat('A',countt);
    accidentalsindividual=resultsmatrix(count,1)*resultsmatrix(count,2)*deltat/timeinterval;
    xlswrite('Tomographyresults.xlsx',[resultsmatrix(count,1:2),resultsmatrix(count,5),accidentalsindividual],Sheet1,xlrange2);
    count = count +1;

end

fclose(scounter);  % close the serial port after the inner loop ends.
pause(2)
%     if stateindexi~=numofstates
%         close all
%     end 

% Writing the average of the measurement results for the state into "Total
% Results" sheet
countt=num2str((stateindexi-1)*2+3);
xl2range2=strcat('A',countt);
accidentalstotal=sum(resultsmatrix(1:numofmeasurements,1))*sum(resultsmatrix(1:numofmeasurements,2))*deltat/(numofmeasurements*timeinterval);
xlswrite('Tomographyresults.xlsx',[sum(resultsmatrix(1:numofmeasurements,1)),sum(resultsmatrix(1:numofmeasurements,2)),sum(resultsmatrix(1:numofmeasurements,5)),accidentalstotal],Sheet2,xl2range2)
pause(statepause);
count=1;

end
clear s;
%% writing date and time of the results into the excle files

timeheader={'year','month','day','hour','minute','seconds'};

% "Gradual Results" sheet
countt=num2str(numofstates*(numofmeasurements+1)+2);
xlrange2=strcat('A',countt);
xlswrite('Tomographyresults.xlsx',timeheader,Sheet1,xlrange2)
countt=num2str(numofstates*(numofmeasurements+1)+3);
xlrange2=strcat('A',countt);
xlswrite('Tomographyresults.xlsx',clockt,Sheet1,xlrange2)

% "Total Results" sheet
xlrangetimeheader=strcat('A',num2str(numofstates*2+2));
xlrangetime=strcat('A',num2str(numofstates*2+3));
xlswrite('Tomographyresults.xlsx',timeheader,Sheet2,xlrangetimeheader)
xlswrite('Tomographyresults.xlsx',clockt,Sheet2,xlrangetime)

% save to "resultsmatrix.txt" file
save('Tomographyresults.txt','resultsmatrix','-ascii')

%% Quantum State Tomography
%Based on "Measurement of Qubits" by James et al.
% Make the column vector named "n" having all 16 measurements per paper
% n = column matrix of counts of the different states enumerated here:
% 1 HH
% 2 HV
% 3 VV
% 4 VH
% Then run the function `r=QSTK(n)', and r will be the 4x4 density matrix.
% The function will do diagnostics on the density matrix that is obtained
% The routines for the diagnostics were taken from Paul Kwiat's posted
% routines in the website
% KG 8/12
% Modified KG 11/15
% Modified KG 7/17
%%
Itxtnam=strcat(datanam,'.txt');
ftxt=fopen(Itxtnam,'wt'); % text file with output
fprintf(ftxt,'%s\n',datanam); % data name

fprintf(ftxt,'%s\n',dateexp); % date
%datanam='James'
datatype='none';
datathick=0;
DelT=50*10^(-9); % coincidence time
Tdwell = timeinterval; %integration time
Tds=num2str(Tdwell);
fprintf(ftxt,' %s %s\n','Tdwell = ',Tds); % data name
fprintf(ftxt,'%s\n','data');
fprintf(ftxt,'%s\n',': coinc     n1       n2  ');
fprintf(ftxt,'%10.0f %10.0f %10.0f\n',n',n1',n2');
n0 = n;
acc = n1 .* n2 .* DelT/Tdwell*0;
n=n-acc;
ibasicdiagnosis = true;
ifulldiagnosis = false;
format compact
%n=n-acc;
%% Gamma matrices (consistent with Appendix A of james et al)
%
% Pauli matrices 0=4
sig{4}=[1 0; 0 1];
sig{1}=[0 1; 1 0];
sig{2}=[0 -1i; 1i 0];
sig{3}=[1 0; 0 -1];
%
Gam{1}=0.5*kron(sig{4},sig{1});
Gam{2}=0.5*kron(sig{4},sig{2});
Gam{3}=0.5*kron(sig{4},sig{3});
Gam{4}=0.5*kron(sig{1},sig{4});
Gam{5}=0.5*kron(sig{1},sig{1});
Gam{6}=0.5*kron(sig{1},sig{2});
Gam{7}=0.5*kron(sig{1},sig{3});
Gam{8}=0.5*kron(sig{2},sig{4});
Gam{9}=0.5*kron(sig{2},sig{1});
Gam{10}=0.5*kron(sig{2},sig{2});
Gam{11}=0.5*kron(sig{2},sig{3});
Gam{12}=0.5*kron(sig{3},sig{4});
Gam{13}=0.5*kron(sig{3},sig{1});
Gam{14}=0.5*kron(sig{3},sig{2});
Gam{15}=0.5*kron(sig{3},sig{3});
Gam{16}=0.5*kron(sig{4},sig{4});
%% Basis vectors (per table I of James et al)
%
PsiH = [1; 0];
PsiV = [0; 1];
PsiD = 2^(-1/2)*[1; 1];
PsiA = 2^(-1/2)*[1; -1];
PsiR = 2^(-1/2)*[1; -1i];
PsiL = 2^(-1/2)*[1; 1i];
Psi{1}=kron(PsiH,PsiH); % HH
Psi{2}=kron(PsiH,PsiV); % HV
Psi{3}=kron(PsiV,PsiV); % VV
Psi{4}=kron(PsiV,PsiH); % VH
Psi{5}=kron(PsiR,PsiH); % RH
Psi{6}=kron(PsiR,PsiV); % RV
Psi{7}=kron(PsiD,PsiV); % DV
Psi{8}=kron(PsiD,PsiH); % DH
Psi{9}=kron(PsiD,PsiR); % DR
Psi{10}=kron(PsiD,PsiD); % DD
Psi{11}=kron(PsiR,PsiD); % RD
Psi{12}=kron(PsiH,PsiD); % HD
Psi{13}=kron(PsiV,PsiD); % VD
Psi{14}=kron(PsiV,PsiL); % VL
Psi{15}=kron(PsiH,PsiL); % HL
Psi{16}=kron(PsiR,PsiL); % RL
Psi;
%%   M-matrices (Appendix B)
B=zeros(16,16);
for iv = 1:16 % Eq. 3.12 of James et al
    for iu = 1:16
        B(iv,iu) = ctranspose(Psi{iv})*Gam{iu}*Psi{iv};
    end
end
B;
Binv = inv(B);
for iu = 1:16 % using corrected Eq. 3.15
    Mtemp = zeros(4,4);
    for iv = 1:16
        Mtemp = Mtemp+Binv(iv,iu)*Gam{iv};
    end
    M{iu} = Mtemp;
end
%% Linear tomography
  norm=n(1)+n(2)+n(3)+n(4);
  s = n./norm; % James Eq. 5.1
  if ibasicdiagnosis 
      disp('s');
      disp(s);
  end
  Lambdanu=(s/norm); % James Eq. 5.15
  rho = zeros (4, 4);
  for j=1:16 % following Eq. 3.20
    rho = rho + n(j) * M{j};
  end
 
  rho = rho/norm;
  if ibasicdiagnosis 
    disp('density linear');
    disp(rho);
  end
  % \sum_\nu s_\nu / (\sum_\nu^4 n_v) 1-4 are the HH, HV, VH, and VV counts
  %% Optimization
  %
  % First step: find the t-parameters from Eqs. 4.4 and 4.6
  Mat11mi1 = [rho(2,2) rho(2,3) rho(2,4);...
      rho(3,2) rho(3,3) rho(3,4);...
      rho(4,2) rho(4,3) rho(4,4)];
  tin(1) = sqrt(det(rho)/det(Mat11mi1));
  Mat1112mi2 = [rho(3,3) rho(3,4); rho(4,3) rho(4,4)];
  tin(2) = sqrt(det(Mat11mi1)/det(Mat1112mi2));
  tin(3) = sqrt(det(Mat1112mi2)/rho(4,4));
  tin(4) = sqrt(rho(4,4));
  Mat12mi1 = [rho(2,1) rho(2,3) rho(2,4);...
      rho(3,1) rho(3,3) rho(3,4);...
      rho(4,1) rho(4,3) rho(4,4)];
  Mat1122mi2 = [rho(3,3) rho(3,4); rho(4,3) rho(4,4)];
  T21 = det(Mat12mi1)/sqrt(det(Mat11mi1)*det(Mat1122mi2));
  tin(5) = real(T21);
  tin(6) = imag(T21);
  Mat1123mi2 = [rho(3,2) rho(3,4); rho(4,2) rho(4,4)];
  T32 = det(Mat1123mi2)/(sqrt(rho(4,4)*det(Mat1122mi2)));
  tin(7) = real(T32);
  tin(8) = imag(T32);
  T43 = rho(4,3)/sqrt(rho(4,4));
  tin(9) = real(T43);
  tin(10) = imag(T43);
  Mat1223mi2 = [rho(3,1) rho(3,4); rho(4,1) rho(4,4)];
  T31 = det(Mat1223mi2)/(sqrt(rho(4,4))*det(Mat1122mi2));
  tin(11) = real(T31);
  tin(12) = imag(T31);
  T42 = rho(4,2)/sqrt(rho(4,4));
  tin(13) = real(T42);
  tin(14) = imag(T42);
  T41 = rho(4,1)/sqrt(rho(4,4));
  tin(15) = real(T41);
  tin(16) = imag(T41);
  %
% input legal density matrix 
  Tmat = [tin(1)         0                  0                0;...
      tin(5)+1i*tin(6)   tin(2)             0                 0;...
      tin(11)+1i*tin(12) tin(7)+1i*tin(8)   tin(3)            0;...
      tin(15)+1i*tin(16) tin(13)+1i*tin(14) tin(9)+1i*tin(10) tin(4)]; % Eq. 4.4
  Rhophysin = ctranspose(Tmat)*Tmat/trace(ctranspose(Tmat)*Tmat); %Eq. 4.5 
  options = optimset('Display', 'off', 'Jacobian', 'off','MaxFunEvals',16000);
iter = 0;
  [t,resnorm,residual,exitflag,output] = lsqnonlin(@Likelihood,tin,[],[],options,n0,norm,acc);
  Tmat = [t(1)       0              0             0;...
      t(5)+1i*t(6)   t(2)           0             0;...
      t(11)+1i*t(12) t(7)+1i*t(8)   t(3)          0;...
      t(15)+1i*t(16) t(13)+1i*t(14) t(9)+1i*t(10) t(4)];
  rhophys = ctranspose(Tmat)*Tmat/trace(ctranspose(Tmat)*Tmat) % Eq. 4.5
  exitflag
  if ifulldiagnosis 
      output
  end
  %output
  rholin = rho;
  rho = rhophys;
  rhor = (rho+conj(rho))/2;

  %% Diagnostics
  %
  if ibasicdiagnosis 
    disp('Tomography Diagnostics');
    disp('The (real) optimized density matrix is');
    disp(rhor);
  end
% Fidelity
  r_HHpVV=0.5*[1,0,0,1;0,0,0,0;0,0,0,0;1,0,0,1];
  r_HHmVV=0.5*[1,0,0,-1;0,0,0,0;0,0,0,0;-1,0,0,1];
  r_HVpVH=0.5*[0,0,0,0;0,1,1,0;0,1,1,0;0,0,0,0];
  r_HVmVH=0.5*[0,0,0,0;0,1,-1,0;0,-1,1,0;0,0,0,0];
  r_HH=[1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0];
  r_VV=[0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1];
  r_HV=[0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0];
  r_VH=[0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 0];
  Plin=trace(rholin);
  rhor2lin=rholin*rholin;
  P2lin=trace(rhor2lin);
  lambdalin=eig(rholin);
  sumofeigenlin=sum(lambdalin);
%
  Pphys=trace(rhophys);
  rhor2phys=rhophys*rhophys;
  P2phys=trace(rhor2phys);
  lambdaphys=eig(rhophys);
  sumofeigenphys=sum(lambdaphys);
  srho1 = r_HHpVV^0.5;
  fid_HHpVVk = trace((srho1*rhophys*srho1)^0.5)^2; % Eq. 29 in Adelpeter et al
%  fid_HHpVV = trace(rho * r_HHpVV);
  if ibasicdiagnosis 
    disp('LINEAR Tr(\rho)');
    disp(Plin);
    disp('LINEAR Tr(\rho^2)');
    disp(P2lin);
    disp('Eigenvalues of LINEAR \rho');
    disp(lambdalin);
    disp('Sum of eigenvalues of LINEAR');
    disp(sumofeigenlin);
    disp('OPTIMIZED Tr(\rho)');
    disp(Pphys);
    disp('OP Tr(\rho^2)');
    disp(P2phys);
    disp('Eigenvalues of OP \rho');
    disp(lambdaphys);
    disp('Sum of eigenvalues');
    disp(sumofeigenphys);
    disp('Fidelity of OP with |H>|H>+|V>|V>');
    disp(fid_HHpVVk);
  end
 fids=num2str(fid_HHpVVk);
  fprintf(ftxt,' %s %s\n','Fidelity with HH+VV = ',fids); % data name

if ifulldiagnosis  
%  disp('Fidelity with |H>|H>+|V>|V>');
  disp('Fidelity with |+1>|H>-|0>|V>');
  fid_HHmVV = trace(rho * r_HHmVV)
  disp('Fidelity with |+1>|V>+|0>|H>');  
  fid_HVpVH = trace(rho * r_HVpVH)
  disp('Fidelity with |+1>|V>-|0>|H>');  
  fid_HVpVH = trace(rho * r_HVpVH)
  disp('Fidelity with |+1>|H>');  
  fid_HH = trace(rho * r_HH)
  disp('Fidelity with |0>|V>');  
  fid_VV = trace(rho * r_VV)
  disp('Fidelity with |+1>|V>');  
  fid_HV = trace(rho * r_HV)
  disp('Fidelity with |0>|H>');  
  fid_VH = trace(rho * r_VH)
end

%% Linear entropy
  d = 4;
  disp('Linear entropy (0 = pure state, 1 = mixed state)')
  lin_e = d * real(1-trace(rho^2)) / (d-1) %Eq. 5.26
  DelP2=0;
  for iv = 1:16 % Eq. 5.28
      Dpds = 0;
      for iu = 1:16
          Dpds = Dpds + trace(M{iu}*M{iv})*s(iu);
      end
      DelP2 = DelP2 + ((8/3)*abs(Dpds))^2*Lambdanu(iv);
  end
  DelP = sqrt(DelP2) 
  DelPs = num2str(DelP);
      
%% Concurrence
%
  Z = [0 0 0 -1; 0 0 1 0; 0 1 0 0; -1 0 0 0];
  R = rho * Z * rho.' * Z;
  % Getting the eigenvalues and eigenvectors of R and R+
  [right,r]=eig(R);
  r;
  Rp=ctranspose(R);
  [rightt,rt]=eig(Rp);
  rt;
%left = inv(right);
  r=real(diag(r));
  right;
  rt=real(diag(rt));
  rightt;
  zetauns = zeros(4,4);
  xiuns = zeros(4,4);
  xi = zeros(4);
  ra = abs(r);
  % sorting the eigenvalues ra and corresponding eigenvectors for both R (zeta) and
  % R+ (xi), which have the same eigenvalues
  for ie1 = 1:4
      ierr = true;
      rsort(ie1,1)=ra(ie1);
      rsort(ie1,2)=ie1;
      rtsort(ie1,1)=rt(ie1);
      rtsort(ie1,2)=ie1;
%       zetauns(:,ie1)=right(:,ie1);
%       for ie2 = 1:4
%           if abs(rt(ie2) - ra(ie1)) < 0.0001 % assume no degeneracy
%               xiuns(:,ie1) = rightt(:,ie2); % eigenvector of R+ for eigenvalue ra
%               ierr = false;
%           end
%       end
%       if ierr == true
%           disp('Could not match eigenvalues of R and R+');
%           disp(' Affects the error in C');
%       end
  end
  rsorted=sortrows(rsort);
  rtsorted=sortrows(rtsort);
  for ie = 1:4
      ra(ie)=rsorted(ie,1);
      zeta(:,ie)=right(:,rsorted(ie,2));
      xi(:,ie)=rightt(:,rtsorted(ie,2));
  end
  ra;
  zeta;
  zeta1=zeta;
  zeta1(:,[2:4])=[];
  zeta2=zeta;
  zeta2(:,[1,3:4])=[];
  zeta3=zeta;
  zeta3(:,[1:2,4])=[];
  zeta4=zeta;
  zeta4(:,[1:3])=[];
  zeta;
  xi1=xi;
  xi1(:,[2:4])=[];
  xi2=xi;
  xi2(:,[1,3:4])=[];
  xi3=xi;
  xi3(:,[1:2,4])=[];
  xi4=xi;
  xi4(:,[1:3])=[];
 
  xi;
  % implement the error according to Eq. 5.35 of James et al
  DelC2=0;
  for iv=1:16
      sqbracket = 0;
      for ia = 1:4
          if ia == 4
              isign = 1;
          else
              isign = -1;
          end
          if ia == 1
              zetacol=zeta1;
              xicol=xi1;
          elseif ia == 2
              zetacol=zeta2;
              xicol=xi2;
          elseif ia == 3
              zetacol=zeta3;
              xicol=xi3;
          else
              zetacol=zeta4;
              xicol=xi4;
          end
          if ra(ia) > 0.001
            for iu = 1:16
                q=M{iu}*Z*M{iv}.'*Z+M{iv}*Z*M{iu}.'*Z;
                bracket=ctranspose(zetacol)*q*s(iu)*xicol;
                sqbracket = sqbracket + isign*0.5*(ra(ia))^(-0.5)*bracket;
            end
          end
      end
      DelC2 = DelC2 + abs(sqbracket)^2*Lambdanu(iv);  
  end
%   for iv = 1:16
      
  tmp = sort(sqrt(r));
  C = real(tmp(4) - tmp(3) - tmp(2) - tmp(1)); %Eq. 5.32
  disp('The Concurrence (1 = Bell state, 0 = product or mixed state)');
  C = max(C, 0)
  DelC=sqrt(DelC2)
  disp('The Tangle'); 
  T=C^2 % 
  DelT=2*C*DelC %Eq. 5.36
  DelTs=num2str(DelT);
%% display bar graph
% clf;
% colormap(jet)
% figure(1)
%  h1=bar3(r_HHmVV);
%  for k1 = 1:length(h1)
%     zdata = h1(k1).ZData;
%     h1(k1).CData = zdata;
%     h1(k1).FaceColor = 'interp';
%  end
% colorbar
%scaled colormap: sjet
sjet = [       0    0.1250    1.0000
         0    0.1875    1.0000
         0    0.2500    1.0000
         0    0.3125    1.0000
         0    0.3750    1.0000
         0    0.4375    1.0000
         0    0.5000    1.0000
         0    0.5625    1.0000
         0    0.6250    1.0000
         0    0.6875    1.0000
         0    0.7500    1.0000
         0    0.8125    1.0000
         0    0.8750    1.0000
         0    0.9375    1.0000
         0    1.0000    1.0000
    0.0625    1.0000    0.9375
    0.1250    1.0000    0.8750
    0.1875    1.0000    0.8125
    0.2500    1.0000    0.7500
    0.3125    1.0000    0.6875
    0.3750    1.0000    0.6250
    0.4375    1.0000    0.5625
    0.5000    1.0000    0.5000
    0.5625    1.0000    0.4375
    0.6250    1.0000    0.3750
    0.6875    1.0000    0.3125
    0.7500    1.0000    0.2500
    0.8125    1.0000    0.1875
    0.8750    1.0000    0.1250
    0.9375    1.0000    0.0625
    1.0000    1.0000         0
    1.0000    0.9375         0
    1.0000    0.8750         0
    1.0000    0.8125         0
    1.0000    0.7500         0
    1.0000    0.6875         0
    1.0000    0.6250         0
    1.0000    0.5625         0
    1.0000    0.5000         0
    1.0000    0.4375         0
    1.0000    0.3750         0
    1.0000    0.3125         0
    1.0000    0.2500         0
    1.0000    0.1875         0
    1.0000    0.1250         0
    1.0000    0.0625         0
    1.0000         0         0
    0.9375         0         0];
h1=figure(2);
colormap(sjet)
% clims=[-0.5 0.5];
  h2=bar3(rhor);
for k = 1:length(h2)
    zdata = h2(k).ZData;
    h2(k).CData = zdata;
    h2(k).FaceColor = 'interp';
end
axis([0.6 4.4 0 5 -0.2 0.5]); 
datathicks=num2str(datathick);
Inam = strcat(datanam,datatype,datathicks);
Ts=num2str(T);
LEs=num2str(lin_e);
Fids=num2str(fid_HHpVVk);
title(['QST',Inam,' um, Tangle=',Ts,' Lin Entr=',LEs,' Fid=',Fids])
print(h1,'-dbmp',Inam);
colorbar
fprintf(ftxt,' %s %s\n','Tangle = ',Ts,' +/- ',DelTs); % tangle
fprintf(ftxt,' %s %s\n','Linear Entropy = ',LEs,' +/- ',DelPs); % linear entropy
