### Measure a RIN spectrum for each current above lasing threshold

#---INPUT 

rin_sweepName = "RIN_current_sweep";    # sweep name to use
rin_resBW=1.00e9;			# Resolution Bandwidth 
rin_minFreq =0.1e9;                     # RIN Spectrum Display Minimum Frequency
rin_maxFreq=30e9;                       # RIN Spectrum Display Maximum Frequency

#---END INPUT

    
##Setting up the sweep

    
#get currents to sweep over
rin_currentIn = li_out.current(li_out.thresholdcurrentindex+1:length(li_out.current));

#setup sweep

if(length(getsweep(rin_sweepName)) >0){ deletesweep(rin_sweepName); }

rin_parameter1Name ="current";
rin_vecParameter1Values = rin_currentIn;
rin_numpts1=length(rin_vecParameter1Values);

addsweep;
setsweep("sweep", "name", rin_sweepName);
setsweep(rin_sweepName, "type", "Values");
setsweep(rin_sweepName, "number of points", rin_numpts1);

#define parameters

rin_para1 = struct;
rin_para1.name=rin_parameter1Name;
rin_para1.parameter = "::Root Element::dcCurrent";
for(rin_i=1:rin_numpts1){
    rin_command = "rin_para1.value_"+num2str(rin_i)+"="+num2str(rin_vecParameter1Values(rin_i))+";";
    eval(rin_command);
}
addsweepparameter(rin_sweepName, rin_para1);

#define results

rin_result = struct;
rin_result.Name="dc power";
rin_result.Result="::Root Element::PWM_1::dc power";
addsweepresult(rin_sweepName,rin_result);
clear(rin_result);

rin_result = struct;
rin_result.Name="noise spectrum";
rin_result.Result="::Root Element::RFSA_1::signal";
addsweepresult(rin_sweepName,rin_result);
clear(rin_result);

rin_result = struct;
rin_result.Name="PIN current";
rin_result.Result="::Root Element::OSC_2::signal";
addsweepresult(rin_sweepName,rin_result);
clear(rin_result);


# run sweep
runsweep(rin_sweepName);
waituntildone;

    
   
##Get results

rin_myr=getsweepresult(rin_sweepName,"dc power");
rin_ePwr=rin_myr.getattribute(rin_myr.getattribute);

rin_myr=getsweepresult(rin_sweepName,"current");
rin_current=rin_myr.getattribute(rin_myr.getattribute);

rin_myr=getsweepresult(rin_sweepName,"noise spectrum");
rin_eSpec=rin_myr.getattribute(rin_myr.getattribute);
rin_eFreq=rin_myr.getparameter("frequency");

mySz=size(rin_eSpec);
rin_simBW=rin_eFreq(mySz(1)); #simulation BW
rin_specRes=rin_eFreq(2)-rin_eFreq(1);
rin_winLen=round(rin_resBW/rin_specRes);
rin_nFreq=round(mySz(1)*rin_maxFreq/rin_simBW)+rin_winLen;
rin_nSweep=mySz(2);
rin_eRIN=rin_eSpec(1:rin_nFreq,1:rin_nSweep);

for(rin_k=1:rin_nSweep){
   rin_eRIN(1:rin_nFreq,rin_k)=rin_eRIN(1:rin_nFreq,rin_k)/rin_ePwr(rin_k);
}

## Plot raw RIN
rin_myLeg=cell(rin_nSweep);
format short; 
for(rin_k=1:rin_nSweep){ 
   rin_myLeg{rin_k} = "cur="+num2str(rin_current(rin_k))+" A";
   }
format long;
plot(rin_eFreq(1:rin_nFreq),10*log10(rin_eRIN),"frequency(Hz)","RIN(dB/Hz)","RIN Spectrum Raw");
setplot("log10x",1);
setplot("x min",rin_minFreq);
setplot("x max",rin_maxFreq);
legend(rin_myLeg);


## Apply low-pass filter 

#define Gaussian window (|G|^2 3dB at ends)
nwin=rin_winLen;
wg=matrix(nwin);
beta=sqrt(2*log(2))/(nwin-1);
    
for (k=1:nwin) { 
     wg(k) = exp(-beta^2*(k-(nwin-1)/2-1)^2); 
     }

#convolve the power spectrum with window to smooth to rbw
rin_eRINs=matrix(rin_nFreq,rin_nSweep);
w=wg;
for(rin_k=1:rin_nSweep){
    yin=rin_eRIN(1:rin_nFreq,rin_k);
    wsum=sum(w^2);
    arrN=size(yin);
    arrM=size(w);
    N=max(arrN);
    M=max(arrM);
    mM=ceil((M-1)/2);
    yout=matrix(N);
    
    for (i=0:N-1){
        nzp=0;
        for (j=0:M-1){  
    
            if( ( i+1-(mM -j)>0) & ( i+1-(mM-j)< N+1) ){
                  yout(i+1,1) = yout(i+1,1) + yin(i+1-(mM -j))*w(j+1)^2;
                  nzp=nzp+1;
            }              
        }
        if(nzp==M) {
           nf=wsum;
           yout(i+1,1)=yout(i+1,1)/wsum;
        }
        else {
           nf = sum(w(max([0,mM- i])+1:min([M-1,N-1+mM -i])+1)^2);
           yout(i+1,1)=yout(i+1,1)/nf;
        }
    }
       rin_eRINs(1:rin_nFreq,rin_k)=yout;
}



##plot smoothed smoothed RIN

plot(rin_eFreq(1:rin_nFreq),10*log10(rin_eRINs),"frequency(Hz)","RIN(dB/Hz)","RIN Spectrum Smoothed");
setplot("log10x",1);
setplot("x min",rin_minFreq);
setplot("x max",rin_maxFreq);
legend(rin_myLeg);