clf
clear

mod=input('What module? (omit the M): ','s');
T=input('What temperature (149,155,or 161)? : ');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load tweaked version of Margie's fringe model parameters
load fringe_data.txt  % Note: this is sorted by fast model ID, not LM ID
fchan_amp  =fringe_data(:,6); % Amplitude (0 to 1)
fchan_phase=fringe_data(:,7); % Phase (cm-1)
fchan_per  =fringe_data(:,8); % Period (cm)
beta       =fringe_data(:,9); % Temp. coef of phase for each EF (cm-1 per K?)
clear fringe_data

% Margie's reference entrance filter temperature
Tefref=154.8020;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load unmodified version of Margie's fringe model parameters
load fringe_data_margie.txt
fchan_amp_margie  =fringe_data_margie(:,6); % Amplitude (0 to 1)
fchan_phase_margie=fringe_data_margie(:,7); % Phase (cm-1)
clear fringe_data_margie


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Info about the modules

% Module names
modules={'12', '11', '10', '9', '8', '7', '6', '5', '4d', '4c', '3', ...
   '4b', '4a', '2b', '1b', '2a', '1a'};

% Number of channels in each module
nchans=[   130, 144,   167, 167, 161, 167, 167, 159,  106,  94,  192, ...
    106, 104,   150,  130,  116,  118];

% Indices (in 1:2378) of each module
indexhi=cumsum(nchans);
indexlo=indexhi - nchans + 1;

% Enterance filter used by each module
filters=[11, 11, 10, 10, 9, 9, 8, 8, 7, 7, 6, 5, 5, 4, 3, 2, 1];
% EF1  = 1a
% EF2  = 2a
% EF3  = 1b
% EF4  = 2b
% EF5  = 4a & 4b
% EF6  = 3
% EF7  = 4c & 4d
% EF8  = 5 & 6
% EF9  = 7 & 8
% EF10 = 9 & 10
% EF11 = 11 & 12


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Freq scaling factor to convert GM freqs (at Tgmref) into freqs
% close to the values in Margie "metrics_155K.txt" file.
%
% Reference grating model temperature
% Tgmref=155.1325;
%
% Note: This gamma is for a different purpose than Margie's gamma!
%
gamma=1./(1 + 3.04E-5);  % eyeballed value

% gamma=2200./2200.065;  % Our GM calibration [~=1/(1 + 2.95E-5)]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load sample SRF measurement data

it=round(T);
if (it == 149)
   Tgm=148.82;
   Tef=148.59;
   eval(['load ../srf/t21621/m' mod '_srf21621.mat'])
elseif (it == 155)
   Tgm=155.15;
   Tef=154.80;
   eval(['load ../srf/t21480/m' mod '_srf21480.mat'])
elseif (it == 161)
   Tgm=160.97;
   Tef=160.45;
   eval(['load ../srf/t21418/m' mod '_srf21418.mat'])
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calc the grating model
[junk, freq, junk, junk, junk, width] = gmodelall2(Tgm, 0);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert module name to index in "modules"
imod=strmatch(mod,modules);
if (length(imod) == 0)
   error ([mod ' is an invalid module; quitting'])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the optimal multiplier of measured SRF freq to match GM

[indgood]=get_good(srf, s.noisearray, s.indok, 2E-2, 9);
[fmid, junk, junk]=measure_width(s.freq, srf, indgood, 0.5);

krev=indexhi(imod):-1:indexlo(imod);  % Pull out channels for module
% Note: SRF data indices are reversed compared to fast model ordering!
fgm=freq(krev)';

junk=fgm(indgood)./fmid(indgood) - 1;
junk2=mean(junk);
junk3=abs( junk - junk2 );
ind=find( junk3 < 2*std(junk) );

fmeasmult2=mean(junk(ind)) + 1;
clear junk2 junk3 ind

% Note: our GM calibration used 2200.065/2200 ~= 1 + 2.95E-5
fmeasmult=1 + 2.95E-5;

plot(indgood,fgm(indgood)./(fmeasmult*fmid(indgood)) - 1,'bo-', ...
     indgood,fgm(indgood)./(fmeasmult2*fmid(indgood)) - 1,'rx-'),grid
title('measured channel centers vs. grating model')

legend({'using our calibration value', 'using optimal value for this data'},0)

pause

% Re-assign using optimal multipler
fmeasmult=fmeasmult2;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calc model SRF; same shape for all channels in module
x=-3:0.01:3; % distance from chan center in units of FWHM
[gfrac, gexp, gslope, lexp, xresid, yresid]=get_srf_params(mod);
[SRF]=zx_srf_resid(x, 0, 1, gfrac, gexp, gslope, lexp, xresid, yresid);

% Loop over channels in this module
kall=indexlo(imod):indexhi(imod);
for k=kall
   ii = k - indexlo(imod) + 1;  % index of channel within module
   iirev=nchans(imod) - ii + 1;

   % Fringe paramters for this channel
   A=fchan_amp(k);      % amplitude
   dch=fchan_per(k);    % period
   pha=fchan_phase(k);  % phase

   % Margie's unmodified fringe paramters for this channel
   Am=fchan_amp_margie(k);      % amplitude
   pham=fchan_phase_margie(k);  % phase

   % Grating Model channel freqs and widths
   f0=freq(k);
   w=width(k);

   % Correction for uncalibrated measured SRF freqs
   fmeas=s.freq(:,iirev)*fmeasmult;
   ymeas=srf(:,iirev);

   % Compute fringes
   v=x*w + f0; % convert x to absolute frequency
   fch=1 - A + A*sin( 2*pi*dch*( v*gamma * ...
      (1 + beta(k)*(Tef - Tefref)) - pha ) ); % calc the fringes
   y=fch.*SRF; % apply the fringes to the model SRF

   % Margie's unmodifed fringes
   fchm=1 - Am + Am*sin( 2*pi*dch*( v*gamma * ...
      (1 + beta(k)*(Tef - Tefref)) - pham ) ); % calc the fringes
   ym=fchm.*SRF; % apply the fringes to the model SRF

   figure(1)  % SRF
   xmeas=(fmeas - f0)./w;
   ind=find(abs(xmeas) < 2.7);
   plot(x,SRF,'k',x,y./max(y),'r',xmeas(ind),ymeas(ind),'c',...
      x,ym./max(ym),'b--'),grid
   title(['M' mod ' index ' int2str(ii)])
   axis([-0.6 0.6 0.5 1.0]);

   legend({'model without fringes', ...
           'model with tweaked fringes', ...
           'measured SRF', ...
           'model with unmodified fringes'}, 0)

   figure(2) % Amplitude and phase
   subplot(2,1,1)
   plot(kall,fchan_phase(kall),'r',k,fchan_phase(k),'ko',...
      kall,fchan_phase_margie(kall),'b--')
   ij=find(fchan_phase(kall) > 0);
   axis([ min(kall), max(kall), min(fchan_phase(kall(ij))), ...
      max( fchan_phase(kall(ij)) )]);
   subplot(2,1,2)
   plot(kall,fchan_amp(kall),'r',k,fchan_amp(k),'ko', ...
      kall,fchan_amp_margie(kall),'b--')
   ij=find(fchan_amp(kall) > 0);
   axis([ min(kall), max(kall), min(fchan_amp(kall(ij))), ...
      max( fchan_amp(kall(ij)) )]);

   legend({'tweaked', 'current channel', 'unmodified'}, 0)

   pause

end % channels

disp(['finished with module ' mod])
