load('Raw_Data.mat') %load raw data; time (ns), voltage traces (mV) and Schottky traces (mW) are kept separate

Time=Time_trace(:,1);

load('LIV_fitted.mat')%LIV.mat is the LIV data fitted with a bi-Gaussian with 5802 data points 
V_b=LIV(:,1);         %bias voltage applied to device (V)
V_d=LIV(:,2);         %voltage measured across the device (mV)
P_meas=LIV(:,3);      %power measured from device by Schottky diode (mW)
                      %corrected for duty cycle (5%) of device operation 
                                        
%approximating the internal field of the device using the voltage across the device

ridge_height=13.9e-6;   %height of QCL ridge
L_p=96e-9;              %spatial-length of the acoustic pulse

%finds the data point where the nominal bias voltage is closest to a value
%in V_b. Then approximates the internal electric field of the device from
%the voltage across the device divided by the ridge height in kV/cm

%pre-allocating space for speed
Voltage=zeros(length(Time),length(bias_voltage));
V_begin=Voltage;
V_smoothed=Voltage;
Electric_Field=Voltage;

r=20;       %value for moving average filter

%approximate internal field at each bias voltage applied to the QCL
[bias_val,bias_index]=min(abs(V_b-bias_voltage)); 
int_field=V_d(bias_index)/ridge_height/100000;

for i=1:length(int_field)
V_begin(:,i)=(Voltage_trace(:,i))/1000/3.3;     
 %convert voltage trace into V and acount for 1x amplifier on measurement
for p=1+r:length(Time)-r-1
    V_smoothed(p,i)=mean(V_begin(p-r:p+r,i));     %smoothing of data
end
Voltage(:,i)=(V_smoothed(:,i)+V_d(bias_index(i)));  %approximated voltage trace in V at each bias voltage
Electric_Field(:,i)=V_smoothed(:,i)/L_p/100000+int_field(i);
%convert voltage trace data in V to electric field across QCL in kV/cm at 
%each bias voltage
end 

%convert measured voltage across device into electric field using thickness
%of QCL ridge
E_d=V_d/13.9e-6/100000;   

% pre-allocating for speed
power_1=zeros(length(Time),length(bias_voltage));
power_offset_removed=power_1;
Power_estimated=power_1;
power_offset_removed_2=power_1;
minValue=power_1;
closestIndex=power_1;

% This for loop finds the points at which the approximated electric field
%(from voltage trace) across the device is closest to a value from E_d
%calculated from the LIV data. It then puts the corresponding power value
%at that index into an array. Offset is removed and expected power value is
%added at each bias voltage.

for x=1:length(bias_voltage)
for j=1:length(Voltage)
    [minValue(j,x),closestIndex(j,x)]=min(abs(E_d-Electric_Field(j,x))); %finds the THz emission power for the
    power_1(j,x)=P_meas(closestIndex(j,x)); % QCL device at that value of internal electric field  
                                            
end
    power_offset_removed(:,x)=power_1(:,x)-mean(power_1(1:1000,x)); %remove any offset from the data
    [P_offset,P_index(x)]=min(abs(E_d-int_field(x))); %Find THz emission power of unperturbed device at each bias voltage
    Power_estimated(:,x)=(power_offset_removed(:,x)*1/88)+P_meas(P_index(x));
    %perturbation is over one period of the structure so divide by no. of periods in whole structure (88)
    %add unperturbed THz emission power (at each bias voltage) to estimated power data
end


%pre-allocating for speed
Schottky=zeros(length(Time),length(bias_voltage));
S_data=Schottky;       
S_smooth=Schottky;

u=20;       %value for moving average filter

for n=1:length(bias_voltage)
    
    S_data(:,n)=Schottky_trace(:,n)*30/10/11.4;  %remove offset
    %Schottky trace measured using AC output, LIV measured using DC output 
    %correction factor of 30 needed. 2x amplifier used so divide trace by 
    %10. Divide by 11.4 for calibration between mV and mW.   
       
for w=1+u:length(Time)-u-1
   S_smooth(w,n)=mean(S_data(w-u:w+u,n));  %smoothing of data 
   
end
   Schottky(:,n)=P_meas(P_index(n)) + S_smooth(:,n);
   %multiply P_meas x 20 to account for duty cycle correction
end
%%

%plot data in nice graphs
for y=1:length(int_field)
    %set time region of interest (ns)
    start(y)=find(Time(:)==20);
    fin(y)=find(Time(:)==120);
    figure(y)
    subplot(2,1,1)
    plot(Time(start(y):fin(y)),Electric_Field(start(y):fin(y),y),'-k')
    xlabel('Time (ns)')
    ylabel('Electric Field (kV/cm)')
    legend('Voltage trace')
    title(['Bias voltage is ' , num2str(bias_voltage(y)) , 'V. Internal field is ' , num2str(round(int_field(y),2)), 'kV/cm'])
    subplot(2,1,2)
    plot(Time(start(y):fin(y)),Power_estimated(start(y):fin(y),y),'-r',Time(start(y):fin(y)),Schottky(start(y):fin(y),y),'-b')
    xlabel('Time (ns)')
    ylabel('Change in Power (mW)')
    legend('Estimated power modulation','Schottky signal')
end

