startup
%SEG=load('DARKSEG1234.txt');
%ASC=load('DARKAS_C.txt');
%POP=load('DARKPOP_LF.txt');
%steprsp=load('step.txt');
stim=steprsp(:,1);
sval=interp1(stim,steprsp(:,2),SEG(:,1));
cr=(1-sval);



dt=1.99968;
dt=dt-2;
t1=0;
t2=30;
tt=ASC(:,1);
ind=find(and(tt>=t1-dt,tt<t2-dt));
ta=ASC(ind,1)+dt;
a=ASC(ind,2)*1000+0.04;



tt=SEG(:,1);
ind=find(and(tt>=t1,tt<t2));
ts=SEG(ind,1);
ssa=interp1(ta,a,ts);

g=-0.06666;
s1=g*SEG(ind,2);
s2=g*SEG(ind,3);
s3=g*SEG(ind,4);
s4=g*SEG(ind,5);

m1=mean(s1(end/4:end))-.5;
m2=mean(s2(end/4:end))-1;
m3=mean(s3(end/4:end))+1.5;
m4=mean(s4(end/4:end));
mm=m1+m2+m3+m4;

s1=s1+54.375*cr *m1/mm;
s2=s2+54.375*cr *m2/mm;
s3=s3+54.375*cr *m3/mm;
s4=s4+54.375*cr *m4/mm;

ss=s1+s2+s3+s4;

t1=.9;
t2=1.5;
ind=find(and(ts>=t1,tt<t2));
ts=ts(ind);
s1=s1(ind);
s2=s2(ind);
s3=s3(ind);
s4=s4(ind);
ss=ss(ind);
ssa=ssa(ind);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find peaks of 00 using the quadrant that sees most 00 (s4)
if 1
    [pk,vpk,pk2,vpk2]=pkFind(ts,s4,0,0,0.8);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find peaks of 10 in ss
[pkss]=pkFind(ts,ss,0,0,0.8);
[pk6]=pkFind(ts,ss,0,0,0.8,600);
pkss=setdiff(pkss,pk6);
%pkss=pkss([1:7,9,11,13,15,17,18]);
%pkss=pkss([5:7,9,11,13,15,17]);
pkss=pkss([7,9,11,13,15,17]);
plot(ts,ss,'-b',pkss,interp1(ts,ss,pkss),'dr',pk,interp1(ts,ss,pk),'dm')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find peaks of 10 in s2
[pks2]=pkFind(ts,s2,0,0,0.3);
%pks2=pks2([2,4:8,10,13,16,19,21,24,26]);
%pks2=pks2([7:8,10,13,16,19,21,24]);
pks2=pks2([10,13,16,19,21,24]);
plot(ts,ss,'-b',pkss,interp1(ts,ss,pkss),'dk',pk,interp1(ts,ss,pk),'dm',ts,s2,'-r',pks2,interp1(ts,s2,pks2),'dk',pk,interp1(ts,s2,pk),'dm')
grid on
xlabel('time (sec)')
ylabel('power (cts)')
legend('AS\_C\_SUM','10','00','AS\_C\_SEG2')
axis([1.2,1.38,-.1,60])
if 0


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fit amplitude shape (not used)
val=spline(pk,interp1(ts,ss,pk),ts);
plot(ts,ss,ts,s1,ts,s2,ts,s3,ts,s4,ts,val);
legend('sum','1','2','3','4','val');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fit the SRM motiopn with a spline and cast in degree
ct=1:length(pk);
cct=spline(pk,ct,ts);
ind=find(and(ts>=1.2,ts<1.38));
deg=360*mod(-cct,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% statistics on 10 peaks
pkssdeg=360*mod(-spline(pk,ct,pkss),1);
pks2deg=360*mod(-spline(pk,ct,pks2),1);
pkssdeg_m=mean(pkssdeg);
pks2deg_m=mean(pks2deg);

pkssdeg_std=std(pkssdeg);
pks2deg_std=std(pks2deg);
pkssdeg_stdm=pkssdeg_std/sqrt(length(pkssdeg));
pks2deg_stdm=pks2deg_std/sqrt(length(pks2deg));

rtgp=mean([pkssdeg_m,pks2deg_m]);
rtgp_state=max(pkssdeg_stdm,pks2deg_stdm);
rtgp_syste=abs(pkssdeg_m-pks2deg_m);

plot(deg(ind),ss(ind),'db',deg(ind),s2(ind),'dr',[pkssdeg_m,pkssdeg_m],[0,55],'-m',[pks2deg_m,pks2deg_m],[0,55],'-m',2*[pkssdeg_m,pkssdeg_m],[0,55],'-m',2*[pks2deg_m,pks2deg_m],[0,55],'-m');
legend('AS\_C\_SUM','AS\_C\_SEG2')
xlabel('Gouy phase (deg, round trip)')
ylabel('power (cts)')
title(['SRC one-way Gouy phase measurement: ',num2str(rtgp/2,4),'deg +-',num2str(rtgp_state/2,2),'deg (statistical error) +-',num2str(rtgp_syste/2,2),'deg (systematic error)']);
grid on
axis([0,180,0,55])


end

