Dear colleague,

   We have written up a paper entitled "Precision formula for the relativistic
corrections to the kinematical Sunyaev-Zeldovich effect for clusters of 
galaxies" and uploaded on astro-ph/0507466.  We hope our results would be of
some use for the upcoming projects for the precision Sunyaev-Zeldovich effect
observations. We would most like to welcome your comments on our paper.

   We have constructed a subroutine that reproduces the results of the present
work.  We are herewith sending it to you for your use.

With my best wishes,


Naoki Itoh
Department of Physics
Sophia University
7-1, Kioi-cho, Chiyoda-ku, Tokyo, 102-8554, Japan
Email: n_itoh@sophia.ac.jp
Telephone: 81 (0)3 3238 3431
Fax:           81 (0)3 3238 3341

********************************************************************
      Implicit Real*8 (A-H, O-Z)

      DMASS = 511.0D0
      Te = 10.D0
      TH = Te/DMASS
      CO = 1.D0
      BETA = 1.D0/300.D0

      DO 10 X = 1.D0/10.D0, 15.D0, 1.D0/10.D0
      SZKIN = SZK(TH,X,CO,BETA)
      SZTH  = SZT(TH,X)
      Write(6, 100) X, SZTH, SZKIN
100   FORMAT(1H , 3(D12.5, 2X))
10    CONTINUE

      STOP
      END
********************************************************************
C===================================================================C
C
C     The function SZT returns the spectral intensity Delta I/y
C     of the thermal SZ effect as a function of (TH, X).
C
C     The function SZK returns the spectral intensity Delta I/y
C     of the kinematic SZ effect as a function of (TH, X, CO, BETA).
C
C     Inputs:
C             TH   = kB Te/mc**2         (real*8)
C             X    = hbar w/kB T         (real*8)
C             BETA = v/c                 (real*8)
C             CO   = cosine(theta_gamma) (real*8)
C
C            Applicable region for TH , X , BETA and CO
C                   0.0 <=  TH  <= 0.05
C                   0.0 <=   X  <= 15.0
C                   0.0 <= BETA <= 1.0/300.0
C                  -1.0 <=  CO  <= 1.0
C     Output:
C              Delta I/y (real*8)
C
C     Remarks: y in this routine is tau (optical depth) in the
C              standard notation
C              v is the peculiar velocity of the cluster
C              c is the velocity of light
C              theta_gamma is the angle between the directions of 
C              the peculiar velocity of cluster and the initial photon
C              momentum (positive z-direction)
C
C                                  Created by Y. Suda (July, 2005)
C
C===================================================================C
      REAL FUNCTION SZT*8(TH, X)
	real*8 X,XX,S2,Y0,Y1,Y2,Y3,Y4,C1,C2,C3,D0,D1,D2,D3,FFF2
	real*8 TH,CO,BETA,FFF1
	
      XX = X*(EXP(+X) + 1.D0)/(EXP(+X) - 1.D0)
      S2 = 4.D0*X**2/(EXP(+X)+EXP(-X)-2.D0)
      
      Y0 =  -4.D0 +            XX
      Y1 = -10.D0 + 47.D0/2.D0*XX - 42.D0/5.D0*XX**2 
     *                            + 7.D0/10.D0*XX**3
     *     + 7.D0/5.D0*S2*(-3.D0 + XX)
      Y2 = -15.D0/2.D0 + 1023.D0/8.D0*XX 
     *     - 868.D0/5.D0*XX**2 + 329.D0/5.D0*XX**3 
     *    -  44.D0/5.D0*XX**4 + 11.D0/30.D0*XX**5
     *     + S2   /30.D0*(-2604.D0 + 3948.D0*XX - 1452.D0*XX**2
     *                             +  143.D0*XX**3)
     *     + S2**2/60.D0*( -528.D0 +  187.D0*XX)
      Y3 = +15.D0/2.D0 + 2505.D0/8.D0*XX 
     *                 - 7098.D0/5.D0*XX**2 +  14253.D0/10.D0*XX**3
     *               - 18594.D0/35.D0*XX**4 + 12059.D0/140.D0*XX**5
     *               -   128.D0/21.D0*XX**6 +    16.D0/105.D0*XX**7
     *     + S2*(-7098.D0/10.D0 + 14253.D0/5.D0*XX 
     *           - 102267.D0/35.D0*XX**2 + 156767.D0/140.D0*XX**3
     *           -    1216.D0/7.D0*XX**4 +       64.D0/7.D0*XX**5)
     *      + S2**2*(-18594.D0/35.D0 + 205003.D0/280.D0*XX
     *               - 1920.D0/7.D0*XX**2 +  1024.D0/35.D0*XX**3)
     *      + S2**3*(-544.D0/21.D0 + 992.D0/105.D0*XX)
      Y4 = -135.D0/32.D0 + 30375.D0/128.D0*XX
     *              -  62391.D0/10.D0*XX**2 + 614727.D0/40.D0*XX**3
     *              - 124389.D0/10.D0*XX**4 + 355703.D0/80.D0*XX**5
     *              -  16568.D0/21.D0*XX**6 +  7516.D0/105.D0*XX**7
     *              -      22.D0/7.D0*XX**8 +    11.D0/210.D0*XX**9
     *      + S2**4*(-682.D0/7.D0 + 7601.D0/210.D0*XX)
     *      + S2**3*(-70414.D0/21.D0 + 465992.D0/105.D0*XX
     *               - 11792.D0/7.D0*XX**2 +  19778.D0/105.D0*XX**3)
     *      + S2**2*(-124389.D0/10.D0 + 6046951.D0/160.D0*XX
     *               - 248520.D0/7.D0*XX**2 + 481024.D0/35.D0*XX**3
     *               -  15972.D0/7.D0*XX**4 + 18689.D0/140.D0*XX**5)
     *      + S2   *(-62391.D0/20.D0 + 614727.D0/20.D0*XX
     *            - 1368279.D0/20.D0*XX**2 + 4624139.D0/80.D0*XX**3
     *            -   157396.D0/7.D0*XX**4 +   30064.D0/ 7.D0*XX**5
     *            -     2717.D0/7.D0*XX**6 +   2761.D0/210.D0*XX**7)  

      FFF1 = X**3/(EXP(+X)-1.D0)
	FFF2 = X*EXP(+X)/(EXP(+X) - 1.D0)
      
	SZT=FFF1*(TH * FFF2 * (Y0 + TH*Y1 + TH**2*Y2 + TH**3*Y3 +TH**4*Y4)
     *)    
      END


      REAL FUNCTION SZK*8(TH, X, CO, BETA)
	real*8 X,XX,S2,Y0,Y1,Y2,Y3,Y4,C1,C2,C3,D0,D1,D2,D3,FFF2
	real*8 TH,CO,BETA,FFF1
	
      XX = X*(EXP(+X) + 1.D0)/(EXP(+X) - 1.D0)
      S2 = 4.D0*X**2/(EXP(+X)+EXP(-X)-2.D0)
      
      Y0 =  -4.D0 +            XX
      Y1 = -10.D0 + 47.D0/2.D0*XX - 42.D0/5.D0*XX**2 
     *                            + 7.D0/10.D0*XX**3
     *     + 7.D0/5.D0*S2*(-3.D0 + XX)
      Y2 = -15.D0/2.D0 + 1023.D0/8.D0*XX 
     *     - 868.D0/5.D0*XX**2 + 329.D0/5.D0*XX**3 
     *    -  44.D0/5.D0*XX**4 + 11.D0/30.D0*XX**5
     *     + S2   /30.D0*(-2604.D0 + 3948.D0*XX - 1452.D0*XX**2
     *                             +  143.D0*XX**3)
     *     + S2**2/60.D0*( -528.D0 +  187.D0*XX)
      Y3 = +15.D0/2.D0 + 2505.D0/8.D0*XX 
     *                 - 7098.D0/5.D0*XX**2 +  14253.D0/10.D0*XX**3
     *               - 18594.D0/35.D0*XX**4 + 12059.D0/140.D0*XX**5
     *               -   128.D0/21.D0*XX**6 +    16.D0/105.D0*XX**7
     *     + S2*(-7098.D0/10.D0 + 14253.D0/5.D0*XX 
     *           - 102267.D0/35.D0*XX**2 + 156767.D0/140.D0*XX**3
     *           -    1216.D0/7.D0*XX**4 +       64.D0/7.D0*XX**5)
     *      + S2**2*(-18594.D0/35.D0 + 205003.D0/280.D0*XX
     *               - 1920.D0/7.D0*XX**2 +  1024.D0/35.D0*XX**3)
     *      + S2**3*(-544.D0/21.D0 + 992.D0/105.D0*XX)
      Y4 = -135.D0/32.D0 + 30375.D0/128.D0*XX
     *              -  62391.D0/10.D0*XX**2 + 614727.D0/40.D0*XX**3
     *              - 124389.D0/10.D0*XX**4 + 355703.D0/80.D0*XX**5
     *              -  16568.D0/21.D0*XX**6 +  7516.D0/105.D0*XX**7
     *              -      22.D0/7.D0*XX**8 +    11.D0/210.D0*XX**9
     *      + S2**4*(-682.D0/7.D0 + 7601.D0/210.D0*XX)
     *      + S2**3*(-70414.D0/21.D0 + 465992.D0/105.D0*XX
     *               - 11792.D0/7.D0*XX**2 +  19778.D0/105.D0*XX**3)
     *      + S2**2*(-124389.D0/10.D0 + 6046951.D0/160.D0*XX
     *               - 248520.D0/7.D0*XX**2 + 481024.D0/35.D0*XX**3
     *               -  15972.D0/7.D0*XX**4 + 18689.D0/140.D0*XX**5)
     *      + S2   *(-62391.D0/20.D0 + 614727.D0/20.D0*XX
     *            - 1368279.D0/20.D0*XX**2 + 4624139.D0/80.D0*XX**3
     *            -   157396.D0/7.D0*XX**4 +   30064.D0/ 7.D0*XX**5
     *            -     2717.D0/7.D0*XX**6 +   2761.D0/210.D0*XX**7)

      C1=10.D0 + (7.D0*S2)/10.D0 - (47.D0*XX)/5.D0 + 
     *  (7.D0*XX**2)/5.D0
	C2=25.D0 + (11.D0*S2**2)/10.D0 - (1117.D0*XX)/10.D0 + 
     *  (847.D0*XX**2)/10.D0 - (183.D0*XX**3)/10.D0 + 
     *  (11.D0*XX**4)/10.D0 + S2*(847.D0/20.D0 - (183.D0*XX)/5.D0 + 
     *  (121.D0*XX**2)/20.D0)
	C3=75.D0/4.D0  + (272.D0*S2**3)/105.D0 - (21873.D0*XX)/40.D0 + 
     *  (49161.D0*XX**2)/40.D0 - (27519.D0*XX**3)/35.D0 + 
     *  (6684.D0*XX**4)/35.D0 - (3917.D0*XX**5)/210.D0 + 
     *  (64.D0*XX**6)/105.D0 + S2**2*(6684.D0/35.D0 - 
     *  (66589.D0*XX)/420.D0 + (192.D0*XX**2)/7.D0) + 
     *  S2*(49161.D0/80.D0 - (55038.D0*XX)/35.D0 + 
     *  (36762.D0*XX**2)/35.D0 -  (50921.D0*XX**3)/210.D0 + 
     *  (608.D0*XX**4)/35.D0)
	C4=-75.D0/4.D0 + (341.D0*S2**4)/42.D0 - (10443.D0*XX)/8.D0 + 
     *  (359079.D0*XX**2)/40.D0 - (938811.D0*XX**3)/70.D0 + 
     *  (261714.D0*XX**4)/35.D0 - (263259.D0*XX**5)/140.D0 + 
     *  (4772.D0*XX**6)/21.D0 - (1336.D0*XX**7)/105.D0 + 
     *  (11.D0*XX**8)/42.D0 + S2**3*(20281.D0/21.D0 - 
     *   (82832.D0*XX)/105.D0 +   (2948.D0*XX**2)/21.D0) + 
     *  S2**2*(261714.D0/35.D0 - (4475403.D0*XX)/280.D0 + 
     *     (71580.D0*XX**2)/7.D0 - (85504.D0*XX**3)/35.D0 + 
     *  (1331.D0*XX**4)/7.D0) + S2*(359079.D0/80.D0 - 
     *   (938811.D0*XX)/35.D0 + (1439427.D0*XX**2)/35.D0 - 
     *     (3422367.D0*XX**3)/140.D0 + (45334.D0*XX**4)/7.D0 -  
     *   (5344.D0*XX**5)/7.D0 + (2717.D0*XX**6)/84.D0)
	D0=-2.D0/3.D0 + (11.D0*XX)/30.D0
	D1=-4.D0 + 12.D0*XX - 6.D0*XX**2 + (19.D0*XX**3)/30.D0 + 
     *   S2*(-3.D0 + (19.D0*XX)/15.D0)
	D2=-10.D0 + (542.D0*XX)/5.D0 - (843.D0*XX**2)/5.D0 + 
     *  (10603.D0*XX**3)/140.D0 - (409.D0*XX**4)/35.D0 + 
     *  (23.D0*XX**5)/42.D0 + 
     *  S2**2*(-409.D0/35.D0 + (391.D0*XX)/84.D0) + 
     *  S2*(-843.D0/10.D0 + (10603.D0*XX)/70.D0 - 
     *     (4499.D0*XX**2)/70.D0 + (299.D0*XX**3)/42.D0)
	D3=-15.D0/2.D0 + (4929.D0*XX)/10.D0 - (39777.D0*XX**2)/20.D0 + 
     *  (1199897.D0*XX**3)/560.D0 - (4392.D0*XX**4)/5.D0 + 
     *  (16364.D0*XX**5)/105.D0 - (3764.D0*XX**6)/315.D0 + 
     *  (101.D0*XX**7)/315.D0 + 
     *  S2**3*(-15997.D0/315.D0 + (6262.D0*XX)/315.D0) + 
     *  S2**2*(-4392.D0/5.D0 + (139094.D0*XX)/105.D0 - 
     *     (3764.D0*XX**2)/7.D0 + (6464.D0*XX**3)/105.D0) + 
     *  S2*(-39777.D0/40.D0 + (1199897.D0*XX)/280.D0 - 
     *     (24156.D0*XX**2)/5.D0 + (212732.D0*XX**3)/105.D0 - 
     *     (35758.D0*XX**4)/105.D0 + (404.D0*XX**5)/21.D0)

      FFF1 = X**3/(EXP(+X)-1.D0)
	FFF2 = X*EXP(+X)/(EXP(+X) - 1.D0)
      SZK =FFF1*(
     *     FFF2 * BETA**2 * (1.D0*Y0/3.D0 
     *      + TH * (5.D0*Y0/6.D0 + 2.D0*Y1/3.D0)
     *      + TH**2 * (5.D0*Y0/8.D0 + 3.D0*Y1/2.D0 + Y2)
     *      + TH**3 *(-5.D0*Y0/8.D0 + 5.D0*Y1/4.D0 + 5.D0*Y2/2.D0 
     *                + 4.D0*Y3/3.D0))
     *    + FFF2 * BETA * CO * (1.D0 + TH*C1 + TH**2*C2 + TH**3*C3 
     *                          + TH**4*C4)
     *    + FFF2 * BETA**2 * (3.D0*CO**2-1.D0)/2.D0 * (D0 + TH*D1 
     *                                           + TH**2*D2 + TH**3*D3))
     
      END
