% mockZehnder.m
% AJW, 1/06

% How can one synthesize a signal
% to be fed into a Pockels cell (PC, or EOM)
% which will produce sidebands at f1 and 5*f1,
% but not at 4*f1 or 6*f1?

% simulate the time series:
Rsamp = 100; % sampling rate, Hz
tmax = 100;  % total time, s
dt = 1/Rsamp;
t = 0:dt:tmax;
fmax = Rsamp/2;
df = 1/tmax;
f = df:df:fmax;
flen = length(f);
fa = (-fmax:df:fmax)-df;
f1 = 1;
f2 = 2*f1;
f3 = 3*f1;
f4 = 4*f1;
f5 = 5*f1;
f6 = 6*f1;

% hefty modulation depths
m1 = 0.2;
m5 = 0.2;
A0 = 1;

% two PC in series
As = A0*exp(i*m1*cos(2*pi*f1*t)).*exp(i*m5*cos(2*pi*f5*t));

% two PC in parallel (Mach Zehnder)
Ap = A0*(exp(i*m1*cos(2*pi*f1*t))+exp(i*m5*cos(2*pi*f5*t)))/2;

% "phase" modulation required for two PC in parallel
Sp = -i*log(Ap/A0);

%
% display the ffts in frequency space from -fNyquist to fNyquist:
Asf = fft(As);
Asf = [Asf((flen+1):end) Asf(1:flen)];
Apf = fft(Ap);
Apf = [Apf((flen+1):end) Apf(1:flen)];
Spf = fft(Sp);
Spf = [Spf((flen+1):end) Spf(1:flen)];

%
figure(1)
semilogy(fa,abs(Asf),'b')
xlabel('frequency (Hz)')
ylabel('fft(As)')
title('PCs in series')
axis([-10 10 1e-2 1e4])
grid on

figure(2)
semilogy(fa,abs(Apf),'b')
xlabel('frequency (Hz)')
ylabel('fft(Ap)')
title('PCs in parallel')
grid on
axis([-10 10 1e-2 1e4])

figure(3)
semilogy(fa,abs(Spf),'b')
xlabel('frequency (Hz)')
ylabel('fft(Sp)')
title('Phase signal for PCs in parallel')
grid on
axis([-10 10 1e-2 1e4])

% do some peak finding 
% to get the value of the fft  at the peaks:
if (0)
% get maxima, and values there:
format short g
indf = find(abs(fa)<200);
ffa = fa(indf);

AA = Asf(indf);
aA = abs(AA);
aA1 = circshift(aA,[0 1]);
aA2 = circshift(aA,[0 -1]);
ind = find(aA > aA1 & aA > aA2 );
p1 = ffa(ind)';
p2 = abs(AA(ind)');
p3 = angle(AA(ind)')/(2*pi);
[p1  p2  p3]

AA = Apf(indf);
aA = abs(AA);
aA1 = circshift(aA,[0 1]);
aA2 = circshift(aA,[0 -1]);
ind = find(aA > aA1 & aA > aA2 );
p1 = ffa(ind)';
p2 = abs(AA(ind)');
p3 = angle(AA(ind)')/(2*pi);
[p1  p2  p3]

AA = Spf(indf);
aA = abs(AA);
aA1 = circshift(aA,[0 1]);
aA2 = circshift(aA,[0 -1]);
ind = find(aA > aA1 & aA > aA2 );
p1 = ffa(ind)';
p2 = abs(AA(ind)');
p3 = angle(AA(ind)')/(2*pi);
[p1  p2  p3]
end

% synthesize a modulation waveform using real time series;
% this works!
for I=1:7,   %I=1 is DC, etc
    fn(I) = (I-1)*f1;
    wn(I) = 2*pi*fn(I);
    Mc = cos(wn(I)*t);
    Ms = sin(wn(I)*t);
    % demodulate, get coefficients for modulation waveform:
    cAs(I) = mean(As.*Mc);
    cAp(I) = mean(Ap.*Mc);
    cSp(I) = mean(Sp.*Mc);
    % synthesize a modulation waveform, using cosines only:
    if (I == 1)                % DC
        Sm = cSp(I)*Mc;
    elseif (I == 4)
    else
        Sm = Sm + 2*cSp(I)*Mc; % 1f, 2f, 4f, 5f, 6f
    end
end

% print out complex fourier coefficients 
% for PCs in parallel, and required phase modulation: 
format short g
[fn' 1e4*cAp' 1e4*cSp']

% all Sp coefficients of SIN waves are zero,
% and the coefficients of COS waves
% are purely imaginary (AM) for DC, 2f, 4f, 6f
% and purely real (PM) for 1f and 5f.

% now display the fft of the phase modulation signal
% and the synthesized MZ. It works!
% But you can't do it with pure PM.

Smf = fft(Sm);
Smf = [Smf((flen+1):end) Smf(1:flen)];
%
figure(4)
semilogy(fa,abs(Smf),'b')
xlabel('frequency (Hz)')
ylabel('fft(Sm)')
title('Synthesized phase signal for PCs in parallel')
axis([-10 10 1e-2 1e4])
grid on

Am = A0*exp(i*Sm);
Amf = fft(Am);
Amf = [Amf((flen+1):end) Amf(1:flen)];

figure(5)
semilogy(fa,abs(Amf),'b')
xlabel('frequency (Hz)')
ylabel('fft(Am)')
title('Synthesized MZ')
grid on
axis([-10 10 1e-2 1e4])

% It must be possible to synthesize a signal
% using pure PM, that supresses the 4f and 6f,
% at the cost of generating unbalanced +-nf sidebands.
% we can live with that, because our detuned SEC
% only makes use of one sideband in a pair, anyway.

% Use only sines and cosines. 
% A guess: try using sines in place of AM cosines.
% The following approach doesn't work: 
% -4f and -6f are suppressed, but +4f and +6f are not.

for I=1:7,   %I=1 is DC, etc
    fn(I) = (I-1)*f1;
    wn(I) = 2*pi*fn(I);
    Mc = cos(wn(I)*t);
    Ms = sin(wn(I)*t);
    % demodulate:
    cSp(I) = mean(Sp.*Mc);
    if (I == 1)
        Sm = imag(cSp(I))*Ms;  % DC
    elseif (I == 2 | I == 6)
        Sm = Sm + 2*real(cSp(I))*Mc;  % 1f and 5f
    elseif (I == 5 | I == 7)
        Sm = Sm + 2*imag(cSp(I))*Ms;  % 4f and 6f
    end
end

Smf = fft(Sm);
Smf = [Smf((flen+1):end) Smf(1:flen)];
%
figure(6)
semilogy(fa,abs(Smf),'b')
xlabel('frequency (Hz)')
ylabel('fft(Sm)')
title('Synthesized phase signal attempt, doesnt work')
axis([-10 10 1e-2 1e4])
grid on

Am = A0*exp(i*Sm);
Amf = fft(Am);
Amf = [Amf((flen+1):end) Amf(1:flen)];

figure(7)
semilogy(fa,abs(Amf),'b')
xlabel('frequency (Hz)')
ylabel('fft(Am)')
title('Synthesized MZ attempt, doesnt work')
grid on
axis([-10 10 1e-2 1e4])

% Something like this might be made to work...