format long

A=1/(2^0.5);
a=0.5;
bone=12.0*sqrt(2.0)/25.0;
btwo=7.0*sqrt(2.0)/50.0;
%bone=sqrt(6.0)/4.0;
%btwo=sqrt(2.0)/4.0;
b=sqrt(0.75);

Maxphase=[];
Maxamplitude=[];

F1=[];
FFT1=[];
amplitudeF11=[];
amplitudeF12=[];
phaseF22=[];
F2=[];
FFT2=[];
amplitudeF21=[];
amplitudeF22=[];
phaseF22=[];
%xi=0.10053;
%xi=0.00155895;cannot dephase throughly
%xi=0.0502;cannot dephase throughly
xi=0.00159;%L=10
%xi=0.10053;

alpha=asin(xi)/2.0;
%alpha=asin(1/n^0.5)*360/2/pi/2 angle
z0=a*conj(a)-b*conj(b);

gama=0;
theta=(2.0*pi/10.0/2500.0);

ket00=[1; 0];
ket11=[0; 1];

F11=[];
amplitude2=[];
phase2=[];
phase12501=[];
amplitude100=[];
finalZn1=[];
Sz=0.5.*[1,0; 0,-1];
Sx=0.5.*[0,1; 1,0];
Sy=0.5.*[0,-1i; 1i,0];

spin=[];
Zn=[];
Xn=[];
Yn=[];
p00=[];
p01=[];
p10=[];
p11=[];
pho00=[];
pho01=[];
pho10=[];
pho11=[];
spinup=[];
order=[];
k=[];
an=[];
finalsum1=[];
finalsum=[];
sum=[];
sumn=[];

pho00all=[];
pho01all=[];
pho10all=[];
pho11all=[];
finalpho001=[];
finalpho011=[];
finalpho101=[];
finalpho111=[];
finalpho00sort=[];
Imfinalpho01sort=[];
finalpho10sort=[];
finalpho11sort=[];

cn=[];
finalc1=[];

finala1=[];
finalZn1=[];
finalXn1=[];
finalYn1=[];
finalk1=[];
middlek=[0];
y=0;
j=1;
finalpho0040middle=[];
finalpho00401=[];

beta=1;
n=1;
k=0;
kplus=[];
    
finalpho001=[];
while j<=2000

        pho00=[pho00,0.25];
    %Zn=[Zn,0];
    %an=[an,0.8];
    %pho00=[pho00,0.25];
    %pho01=[pho01,0.-0.25];
    %pho10=[pho10,0.5];
    %pho11=[pho11,0.5];

    while beta<=100000

        %a0=a*cos(gama)-b*1i*sin(gama);
        %b0=b*cos(gama)-a*1i*sin(gama);

        %P=sqrt(a0*conj(a0)+b0*conj(b0));

        %a=a0/P;
        %b=b0/P;
    
        while n<=10

            %driving;
            a0=a*cos(theta)-b*1i*sin(theta);
            b0=b*cos(theta)-a*1i*sin(theta);

            %a=a0
            %b=b0

            P=sqrt(a0*conj(a0)+b0*conj(b0));

            a=a0/P;
            b=b0/P;

            %measure
            %the probability of measuring up
            c0=a*(cos(alpha)+sin(alpha))*A;
            d0=b*(cos(alpha)-sin(alpha))*A;
            %normalize=sqrt(c0*conj(c0)+d0*conj(d0));
            %c=c0/normalize;
            %d=d0/normalize;

            pup=c0*conj(c0)+d0*conj(d0);
    
            s=0+(1-0)*rand;
            if s<=pup;

                a0=a*(cos(alpha)+sin(alpha))*A;
                b0=b*(cos(alpha)-sin(alpha))*A;

                P=sqrt(a0*conj(a0)+b0*conj(b0));

                a=a0/P;
                b=b0/P;
                %when the test spin spins up, the state of the qubit after normalization:
                %phiup=a*ket00+b*ket11;
                %the partial magnetization in Z
                %Z=2*phiup.'*Sz*phiup;
                %Z=(conj(a)*a-conj(b)*b);
                %the partial magnetization in X
                %X=(conj(a)*b+conj(b)*a);
                %the partial magnetization in Y
                %Y=(conj(b)*a-conj(a)*b)*1i;
                %sum=X*conj(X)+Y*conj(Y)+Z*conj(Z);
                %when the test spin spins up, the density matrix of the qubit
                %phoup=phiup.*phiup';
                %components of the density matrix
                %pho00every=a*conj(a);
                %pho01every=a*conj(b);
                %pho10every=phoup(2,1);
                %pho11every=phoup(2,2);

                %spin=[spin,"+1"];
                %Zn=[Zn,Z];
                %y=y+1;
                %k=[k,y];
                %k=k+1;

            
            else

                a0=a*(cos(alpha)-sin(alpha))*A;
                b0=b*(cos(alpha)+sin(alpha))*A;

                P=sqrt(a0*conj(a0)+b0*conj(b0));

                a=a0/P;
                b=b0/P;
 
                %when the test spin spins down, the state of the qubit after normalization:
                %phidown=a*ket00+b*ket11;
                %the partial magnetization in Z
                %Z=conj(a)*a-conj(b)*b;
                %the partial magnetization in X
                %X=conj(a)*b+conj(b)*a;
                %the partial magnetization in Y
                %Y=(conj(b)*a-conj(a)*b)*1i;
                %sum=X*conj(X)+Y*conj(Y)+Z*conj(Z);
                %when the test spin spins down, the density matrix of the qubit
                %phodown=phidown.*phidown';
                %components of the density matrix
                %pho00every=a*conj(a);
                %pho01every=a*conj(b);
                %pho10every=phodown(2,1);
                %pho11every=phodown(2,2);

                %spin=[spin,"-1"];
                %Zn=[Zn,Z];
                %k=[k,0];
                %k=k+0;

            end 

        n=n+1;

        end

    n=1;
    y=0;
    beta=beta+1;
    pho001=a*conj(a);
    pho00=[pho00,pho001];
    %Z=conj(a)*a-conj(b)*b;
    %finalZn1=[finalZn1;Z];
    %kplus=[kplus;k];
    %k=0;
    %Xn=[Xn,X];
    %Yn=[Yn,Y];
    %sumn=[sumn,sum];
    %pho00=[pho00,pho00every];
    %pho01=[pho01,pho01every];
    %pho10=[pho10,pho10every];
    %pho11=[pho11,pho11every];
    
    %an=[an,a];
 
    %temp=0;
    %for idx=1:length(k)
    %    if k(idx)>=temp
    %        temp=k(idx);
    %    end
    %end
    %middlek=[middlek,temp];

    %k=[];

    end

    j=j+1;
    %finala1=[finala1;an];
    %an=[];
    %finalZn1=[finalZn1;Zn];
    %Zn=[];
    %finalXn1=[finalXn1;Xn];
    %Xn=[];
    %finalYn1=[finalYn1;Yn];
    %Yn=[];
    %finalsum1=[finalsum1;sumn];
    %sumn=[];
    %finalk1=[finalk1;middlek];
    %middlek=[0];
    finalpho001=[finalpho001;pho00];
    pho00=[];
    %finalpho011=[finalpho011;pho01];
    %pho01=[];
    %finalpho101=[finalpho101;pho10];
    %pho10=[];
    %finalpho111=[finalpho111;pho11];
    %pho11=[];

     beta=1;
     n=1;
     a=0.5;
     b=sqrt(0.75);
     z0=a*conj(a)-b*conj(b);

end

%finala=finala1.';
%finalasort=sort(finala,'ascend');
%finalk1i=real(finalk1);
%finalk=finalk1i.';
%finalksort=sort(finalk,'ascend');
%signal=finalk./100;
%finalZn=finalZn1.';
finalpho00=finalpho001.';
row_averages = mean(finalpho00, 2);

% Define the segments
%segments = [
%    1, 500;
%    49751, 50250;
%    99501, 100000
%];

% Initialize the output variables
%amplitudes = zeros(4000, 3);
%phases = zeros(4000, 3);

% Iterate over segments
%for s = 1:size(segments, 1)
    % Perform FFT, calculate amplitude and phase, and find maximum amplitude and corresponding phase
%    for i = 1:4000
%        fft_row = fft(finalpho001(i, segments(s, 1):segments(s, 2)));
%        amplitude = abs(fft_row);
%        phase = angle(fft_row);
%
%        [max_amplitude, max_idx] = max(amplitude);
%        amplitudes(i, s) = max_amplitude;
%        phases(i, s) = phase(max_idx);
%    end
%end

% Round amplitudes and phases to integers
%int_amplitudes = round(amplitudes);
%int_phases = round(phases);

% Count the occurrences of different (amplitude, phase) pairs for each segment
%for s = 1:size(segments, 1)
%    unique_pairs = unique([int_amplitudes(:, s), int_phases(:, s)], 'rows');
%    pair_count = size(unique_pairs, 1);
%    %fprintf('There are %d unique (amplitude, phase) pairs in segment %d.\n', pair_count, s);
%end

%xlswrite('heatmap,b=0.6124+0.3536i,N=2500.xlsx',pair_count);
%xlswrite('pho00,a=0.5,N=2500.xlsx',finalpho00);
xlswrite('averagepho00,a=0.5,N=2500.xlsx',row_averages);

%xlswrite('k,a=0.707,b=0.707,N=25.xlsx',kplus);
%finalXn=finalXn1.';
%finalYn=finalYn1.';
%finalsum=finalsum1.';
%finalZnsort=sort(finalZn,'ascend');
%finalpho009=finalpho001.';
%F=F2.';
%amplitude9=amplitude2.';
%amplitude8=amplitude9.*0.001598;
%phase9=phase2.';

%finalpho00=round(finalpho009,2);
%amplitude=round(amplitude8,2);
%phase=round(phase9,2);

%numberamplitude=tabulate(Maxamplitude);
%numberphase=tabulate(Maxphase);

%finalpho00sort=sort(finalpho00,'ascend');
%finalpho01=finalpho011.';
%Imfinalpho01=imag(finalpho01);
%Imfinalpho01sort=sort(Imfinalpho01,'ascend');
%finalpho10=finalpho101.';
%finalpho10sort=sort(finalpho10,'ascend');
%finalpho11=finalpho111.';
%finalpho11sort=sort(finalpho11,'ascend');

%aa=reshape(finala,[],1);
%kk=reshape(finalk,[],1);
%ZnZn=reshape(finalZn,[],1);
%0=reshape(finalpho00,[],1);
%pho01pho01=reshape(Imfinalpho01,[],1);
%pho10pho10=reshape(finalpho10,[],1);
%pho11pho11=reshape(finalpho11,[],1);
%signalsignal=reshape(signal,[],1);