%get and plot responsivity
%responsivity is calculated from derivative to avoid having to know dark
%current
clear all

%Source meter settings:
compliance = 100e-3;    %amps
Vstart = 0;             %V
Vstop = -6;
Vstep = -0.25;
pse = 0.08;             %sec (pause between data points)
Vvec = Vstart:Vstep:Vstop;
Vvec_1 = Vstop:-Vstep:Vstart;

%Attenuator settings:
Pstart = 23;           %dBm
Pstop = 20;
Pstep = -.5;
Apause = 2;             %sec
Pvec = Pstart:Pstep:Pstop;
Loss = 4;               %dB - attenuator-to-PD loss, for respo calculation

%GPIB:
meterGPIB = 30;
attGPIB = 20;

%initialize source meter + set voltage to Vstart
address = strcat('GPIB0::', num2str(meterGPIB), '::INSTR');
source = visa('ni', address);
fopen(source)
voltage = Vvec(1);
fprintf(source, ':SENS:CURR:PROT %e', compliance);
fprintf(source, ':SOUR:VOLT:LEV %e', voltage);
fprintf(source, ':OUTP ON');
fprintf(source, ':READ?')
current_holder = str2num(fscanf(source));
current = current_holder(2);

%initialize attenuator + set power to Pstart:
address2 = strcat('GPIB0::', num2str(attGPIB), '::INSTR');
att = visa('ni', address2);
fopen(att)
fprintf(att, 'OUTP1:CHAN1:POW %d dBm', Pvec(1));
fprintf(att, 'OUTP1 ON')

%Outer loop: optical power:
%current matrix has
    %row = power level;
    %column = set voltage
for j = 1:length(Pvec)
    %set optical power:
    fprintf(att, 'OUTP1:CHAN1:POW %d dBm', Pvec(j))
    pause(Apause);
    %get power in dBm:
    fprintf(att, ':READ1:CHAN1:POW?')
    power(j) = str2num(fscanf(att));
    %step through voltages: go up:
    for i = 1:length(Vvec)
        voltage = Vvec(i);
        fprintf(source, ':SOUR:VOLT:LEV %e; *WAI', voltage);
        fprintf(source, ':READ?; *WAI')
        current_holder = str2num(fscanf(source));
        current(j, i) = current_holder(2);
        pause(pse);
    end
    %go back down slowly
    for i = 1:length(Vvec_1)
        voltage = Vvec_1(i);
        fprintf(source, ':SOUR:VOLT:LEV %e; *WAI', voltage);
        pause(pse);
    end
end

fclose(source)
fclose(att)

uisave()

power_mW = 10.^((power-Loss)/10);
V_plot = -1:-1:ceil(min(Vvec));
for k = 1:length(V_plot)
    [a b] = min(abs(Vvec-V_plot(k)));
    index(k) = b;
end
figure(1)
red = linspace(0, 1, length(V_plot));
blue = linspace(1, 0, length(V_plot));
green = zeros(1, length(V_plot));
plotlist = cat(1, red, green, blue);
for k = 1:length(V_plot)
    H(k)=plot(power_mW, -1*current(:,index(k)));
    set(H(k), 'Color', plotlist(:,k)');
    hold on
end
xlabel('Optical power (mW)')
ylabel('Current (A)')
hold off
figure(2)
for k = 1:length(V_plot)
    Respo = -1e3*diff(current(:,index(k)))./diff(power_mW)';
    H(k)=plot(power_mW(2:length(power_mW)), Respo);
    set(H(k), 'Color', plotlist(:,k)');
    hold on
end
xlabel('Optical power (mW)')
ylabel('Responsivity (A/W)')
hold off

    

