utine harm24(x,a,b)

        implicit none

        real    x(24),a(13),b(12)

c ***** THIS ROUTINE PERFORMS A 24 POINT HARMONIC ANALYSIS ON THE
c ***** VECTOR x. RETURNED ARE THE FOURIER COEFFICIENTS a (COSINE
c ***** TERMS) AND b (SINE TERMS).

        real    u(12),v(12),r(6),s(6),p(6),q(6),zl(3),zm(3),g(3)
        real    h(3),e,f,c,d,cc,ss,aa,bb,v1,v7,ex(24),px(24)
        integer i,j
        equivalence     (px(1),u(1),zl(1))
        equivalence     (px(4),zm(1))
        equivalence     (px(7),g(1))
        equivalence     (px(10),h(1))
        equivalence     (v(2),c)
        equivalence     (v(3),d)
        equivalence     (v(4),cc)
        equivalence     (v(5),ss)
        equivalence     (v(6),aa)
        equivalence     (v(7),bb)
        equivalence     (px(13),v(1))
        equivalence     (ex(1),r(1))
        equivalence     (ex(7),s(1))
        equivalence     (ex(13),p(1))
        equivalence     (ex(19),q(1))

        do i = 2,12
           j    = 26-i
           u(i) = x(i)+x(j)
           v(i) = x(i)-x(j)
        enddo
        u(1) = x(1)+x(13)
        v(1) = x(1)-x(13)

        do i = 2,6
           j    = 14-i
           r(i) = u(i)+u(j)
           s(i) = u(i)-u(j)
        enddo
        r(1) = u(1)+u(7)
        s(1) = u(1)-u(7)

        do i = 2,6
           j    = 14-i
           p(i) = v(i)+v(j)
           q(i) = v(i)-v(j)
        enddo
        v1 = v(1)
        v7 = v(7)

        zl(1) = r(1)+r(4)
        zl(2) = r(2)+r(6)
        zl(3) = r(3)+r(5)
        zm(1) = r(1)-r(4)
        zm(2) = r(2)-r(6)
        zm(3) = r(3)-r(5)
        g(2)  = q(2)+q(6)
        h(2)  = q(2)-q(6)
        g(3)  = q(3)+q(5)
        h(3)  = q(3)-q(5)
        e     = zl(2)+zl(3)
        f     = zl(2)-zl(3)
        c     = h(2)+h(3)
        d     = h(2)-h(3)
        ss = 0.258819e0
        cc = 0.965926e0
        bb = 0.707107e0
        aa = 0.866025e0

c ***** COSINE COEFFS.
        a(1)  = zl(1)+e
        a(2)  = v1+cc*s(2)+aa*s(3)+bb*s(4)+0.5*s(5)+ss*s(6)
        a(3)  = s(1)+aa*zm(2)+0.5*zm(3)
        a(4)  = v1+bb*(s(2)-s(4)-s(6))-s(5)
        a(5)  = zm(1)+0.5*f
        a(6)  = v1+ss*s(2)-aa*s(3)-bb*s(4)+0.5*s(5)+cc*s(6)
        a(7)  = s(1)-zm(3)
        a(8)  = v1-ss*s(2)-aa*s(3)+bb*s(4)+0.5*s(5)-cc*s(6)
        a(9)  = zl(1)-0.5*e
        a(10) = v1-bb*(s(2)-s(4)-s(6))-s(5)
        a(11) = s(1)-aa*zm(2)+0.5*zm(3)
        a(12) = v1-cc*s(2)+aa*s(3)-bb*s(4)+0.5*s(5)-ss*s(6)
        a(13) = zm(1)-f

c ***** SINE COEFFS.
        b(2)  = ss*p(2)+0.5*p(3)+bb*p(4)+aa*p(5)+cc*p(6)+v7
        b(3)  = 0.5*g(2)+aa*g(3)+q(4)
        b(4)  = p(3)-v7+bb*(p(2)+p(4)-p(6))
        b(5)  = aa*c
        b(6)  = cc*p(2)+0.5*p(3)-bb*p(4)-aa*p(5)+ss*p(6)+v7
        b(7)  = g(2)-q(4)
        b(8)  = cc*p(2)-0.5*p(3)-bb*p(4)+aa*p(5)+ss*p(6)-v7
        b(9)  = aa*d
        b(10) = v7-p(3)+bb*(p(2)+p(4)-p(6))
        b(11) = 0.5*g(2)-aa*g(3)+q(4)
        b(12) = ss*p(2)-0.5*p(3)+bb*p(4)-aa*p(5)+cc*p(6)-v7

        do i = 1,12
           a(i) = a(i)/12.e0
           b(i) = b(i)/12.e0
        enddo
        a(13) = a(13)/12.e0

        return
        end