sz_relcor

 f = sz_relcor(freq, kTe)

 this is the sz f(x) frequency dependence with the relativistic
 correction to 5th order. this is what Sam LaRoque used, transcribed
 to matlab.

 freq is in GHz.
 kTe is in keV.

 Itoh el al. 
 ApJ 502:7-15,1998 July 20
 equations 2.8 , 2.26-30, 2.32-33

 http://www.journals.uchicago.edu/ApJ/journal/issues/ApJ/v502n1/36888/36888.html

0001 function f = sz_relcor(freq, kTe)
0002 
0003 % f = sz_relcor(freq, kTe)
0004 %
0005 % this is the sz f(x) frequency dependence with the relativistic
0006 % correction to 5th order. this is what Sam LaRoque used, transcribed
0007 % to matlab.
0008 %
0009 % freq is in GHz.
0010 % kTe is in keV.
0011 %
0012 % Itoh el al.
0013 % ApJ 502:7-15,1998 July 20
0014 % equations 2.8 , 2.26-30, 2.32-33
0015 %
0016 % http://www.journals.uchicago.edu/ApJ/journal/issues/ApJ/v502n1/36888/36888.html
0017  
0018 
0019 if nargin<2
0020   kTe=0;
0021 end
0022 
0023 freq = freq*1e9;
0024 
0025 hP = 6.6260693e-34; % planck in J*s
0026 k = 1.3806505e-23; % J/K = kg m^2 s^-2 / K
0027 me = 9.1093826e-31; % kg
0028 Tcmb = 2.726; % K
0029 c = 299729458; % m/s
0030 keV = 1.1605e7; % K/keV
0031 
0032 x = hP .* freq ./ (k.*Tcmb);
0033 
0034 xb1 = x ./ tanh(x/2); 
0035 % xb1 = x*(exp(x)+1)/(exp(x)-1);
0036 
0037 xb2 = xb1.^2;
0038 xb3 = xb1.^3;
0039 xb4 = xb1.^4;
0040 xb5 = xb1.^5;
0041 xb6 = xb1.^6;
0042 xb7 = xb1.^7;
0043 xb8 = xb1.^8;
0044 xb9 = xb1.^9;
0045 
0046 sb1 = x ./ sinh(x/2);
0047 sb2 = sb1.^2;
0048 sb4 = sb1.^4;
0049 sb6 = sb1.^6;
0050 sb8 = sb1.^8;
0051 
0052 theta1 = (kTe*keV*k)./(me*c^2);
0053 theta2 = theta1.^2;
0054 theta3 = theta1.^3;
0055 theta4 = theta1.^4;
0056 
0057 y0 = xb1 - 4;
0058  
0059 y1 = -10 + (47/2)*xb1 - (42/5)*xb2 + (7/10)*xb3 + sb2.*(-21/5 + (7/5)*xb1);
0060 
0061 y2 = -(15/2) + (1023/8)*xb1 - (868/5)*xb2 + (329/5)*xb3 - (44/5)*xb4 + (11/30)*xb5 ...
0062  + sb2.*(-434/5 + (658/5)*xb1 - (242/5)*xb2 + (143/30)*xb3) + sb4.*(-44/5+(187/60)*xb1);
0063 
0064 y3 = 15/2 + (2505/8)*xb1 - (7098/5)*xb2 + (14253/10)*xb3 - (18594/35)*xb4 + (12059/140)*xb5 - (128/21)*xb6 + (16/105)*xb7 ...
0065  + sb2.*(-7098/10 + (14253/5)*xb1 - (102267/35)*xb2 + (156767/140)*xb3 - (1216/7)*xb4 + (64/7)*xb5)...
0066  + sb4.*(-18594/35 + (205003/280)*xb1 - (1920/7)*xb2 + (1024/35)*xb3) + sb6.*(-544/21 + (922/105)*xb1) ;
0067 
0068 y4 = (-135/32) + (30375/128)*xb1 - (62391/10)*xb2 + (614727/40)*xb3 - (124389/10)*xb4 ...
0069  +(355703/80)*xb5 - (16568/21)*xb6 + (7516/105)*xb7 - (22/7)*xb8 + (11/210)*xb9 ...
0070  + sb2.*( -62391/20 + (614727/20)*xb1 - (1368279/20)*xb2 + (4624139/80)*xb3 -(157396/7)*xb4 ...
0071  + (30064/7)*xb5  - (2717/7)*xb6 + (2761/210)*xb7 ) ...
0072  + sb4.*(-124389/10 + (6046951/160)*xb1 - (248520/7)*xb2 +(481024/35)*xb3- (15972/7)*xb4 + (18689/140)*xb5 ) ...
0073  + sb6.*(-70414/21 + (465992/105)*xb1 - (11792/7)*xb2 + (19778/105)*xb3) ...
0074  + sb8.*(-682/7 + (7601/210)*xb1 );
0075 
0076 f = y0'*ones(size(theta1)) + (y1'*theta1) + (y2'*theta2) + (y3'*theta3) + (y4'*theta4);