function [TDI2Out, TDI1Out] = TDIComb(fs, s21, s31, s12, s13, tau33, tau13, beta3o, tau22, tau12, beta2o)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [TDI2Out, TDI1Out] = TDIComb(fs, s21, s31, s12, s13, tau33,
%       tau13, beta3, tau22, tau12, beta2)
%
% Computes the TDI 1.0 and TDI 2.0 Combinations based on the 6-variable
% arm-length delays. Based on "Data combinations accounting for LISA
% spacecraft motion" by Shaddock et. al.
%
% INPUTS:
%  fs = Sample frequency
%
%  s21 = Sensor signal (sending: SC2 receiving: SC1)
%  s31 = Sensor signal (sending: SC3 receiving: SC1)
%  s12 = Sensor signal (sending: SC1 receiving: SC2)
%  s13 = Sensor signal (sending: SC1 receiving: SC3)
%
%  tau33 = Round-trip delay from SC1 through SC2 (Arm3)
%  tau13 = One-way delay time from SC1 to SC2 (Arm3)
%  beta3 = Arm3 time compression factor
%  beta = (velocity between SC1 & SC2/speed of light)
%  tau22 = Round-trip delay from SC1 through SC3 (Arm2)
%  tau12 = One-way delay time from SC1 to SC3 (Arm2)
%  beta2 = Arm2 time compression factor
%
% OUTPUTS:
%  TDI2Out = TDI2.0 Combination
%  TDI1Out = TDI1.0 Combination
%
%   Author: Shawn Mitryk
%   Created: 9-18-11
%   Updated: 3-28-12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	L1 = length(s12);
	N = 51;
	M = (N-1)/2;

	alpha2 = (1-2*beta2o);
	alpha3 = (1-2*beta3o);
	gamma2 = (1-beta2o);
	gamma3 = (1-beta3o);
	beta3 = 1-alpha3;
	beta2 = 1-alpha2;


% TDI 2.0 shift values
	Tau23 = (tau22)*alpha3+(tau33)
	beta23 = 1-(alpha2*alpha3)
	Tau32 = (tau33)*alpha2+(tau22)
	beta32 = 1-(alpha3*alpha2)
	Tau232 = alpha2^2*tau33 + alpha2*tau22 + tau22
	beta232 = 1-(alpha2^2*alpha3)
	Tau323 = alpha3^2*tau22 + alpha3*tau33 + tau33
	beta323 = 1-(alpha3^2*alpha2)


% Calculate the sensor signals from the one-way signals
	[chan12.data, chan12.t] = TVFDfilter(s12,fs,beta2o,tau12,N);
	[chan13.data, chan13.t] = TVFDfilter(s13,fs,beta3o,tau13,N);
	chan21.data = s21(2*N:L1-M);
	chan31.data = s31(2*N:L1-M);
	chan22.data = chan21.data + chan12.data;
	chan33.data = chan31.data + chan13.data;


% Form the TDI Combinations based on the 6-variable delays
L2 = length(chan22.data);

	[chan2_2s.data,chan2_2s.t] = TVFDfilter(chan22.data,fs,beta2,tau22,N);
	[chan3_3s.data,chan3_3s.t] = TVFDfilter(chan33.data,fs,beta3,tau33,N);
	[chan22_23s.data,chan2_23s.t] = TVFDfilter(chan22.data,fs,beta23,Tau23,N);
	[chan33_32s.data,chan3_32s.t] = TVFDfilter(chan33.data,fs,beta32,Tau32,N);
	[chan22_232s.data,chan2_232s.t] = TVFDfilter(chan22.data,fs,beta232,Tau232,N);
	[chan33_323s.data,chan3_323s.t] = TVFDfilter(chan33.data,fs,beta323,Tau323,N);


	chan22a.data = chan22.data(2*N:L2-M)';
	chan33a.data = chan33.data(2*N:L2-M)';
	chan22_2a.data = chan22_2s.data;
	chan33_3a.data = chan33_3s.data;
	chan22_23a.data = chan22_23s.data;
	chan33_32a.data = chan33_32s.data;
	chan22_232a.data = chan22_232s.data;
	chan33_323a.data = chan33_323s.data;

	tdi1a.data = chan22a.data + chan33a.data - chan22_2a.data - chan33_3a.data;
	tdi1b.data = - chan22_23a.data - chan33_32a.data + chan22_232a.data + chan33_323a.data;
	tdi2.data = tdi1a.data + tdi1b.data;

	TDI1Out = tdi1a.data;
	TDI2Out = tdi2.data;