! File: fft13t.i90
! Public domain 2004 James Van Buskirk
subroutine fft13t(h)
! 164 adds, 64 muls
   complex(wp), intent(inout) :: h(0:12)
   real(wp) xr1, xi1, xr2, xi2, xr3, xi3, xr4, xi4, xr5, xi5
   real(wp) xr6, xi6, xr7, xi7
   real(wp) yr1, yi1, yr2, yi2, yr3, yi3, yr4, yi4, yr5, yi5
   real(wp) yr6, yi6
   real(wp) pr2, pi2, pr3, pi3, pr4, pi4
   real(wp) pr5, pi5, pr6, pi6, pr7, pi7
   real(wp) pr234, pi234, sr1, si1, ur234, ui234
   real(wp) ur2, ui2, ur3, ui3, ur4, ui4
   real(wp) ur5, ui5, ur6, ui6, ur7, ui7
   real(wp) sr2, si2, sr3, si3, sr4, si4
   real(wp) sr5, si5, sr6, si6, sr7, si7
   real(wp) qr1, qi1, qr3, qi3, qr_5, qi_5
   real(wp) qr2, qi2, qr4, qi4, qr_6, qi_6
   real(wp) vr1, vi1, vr2, vi2
   real(wp) qr34, qi34, qr_56, qi_56
   real(wp) vr34, vi34, vr3, vi3, vr4, vi4
   real(wp) vr56, vi56, vr5, vi5, vr6, vi6
   real(wp) tr1, ti1, tr3, ti3, tr5, ti5
   real(wp) tr2, ti2, tr4, ti4, tr6, ti6
   real(wp), parameter :: dc1 = -0.66277335223429382656331008444690382_wp
   real(wp), parameter :: dc2 = 0.73124599097534822519618254560377760_wp
   real(wp), parameter :: dc3 = 1.0070740657275332544937477077369340_wp
   real(wp), parameter :: dc4 = -0.30816846519175820059367219184688543_wp
   real(wp), parameter :: dc5 = 0.81698338691215549726750306085822509_wp
   real(wp), parameter :: dc6 = 0.22376403323791637458136993583748652_wp
   real(wp), parameter :: ds1 = 0.57514072947400312136838554745545335_wp
   real(wp), parameter :: ds2 = -0.17413860115213590500566079492926474_wp
   real(wp), parameter :: ds3 = -0.33582506518644535421963182119524142_wp
   real(wp), parameter :: ds4 = 4.5240494294812713569277280991401412E-0002_wp
   real(wp), parameter :: ds5 = 1.1543953381323634420147226757584967_wp
   real(wp), parameter :: ds6 = -8.7981928766792081008399945211312915E-0002_wp
   real(wp), parameter :: ds7 = 0.90655220171271016880349079977456841_wp
   real(wp), parameter :: ds8 = 1.1971367726043428094538453399784083_wp
   real(wp), parameter :: ds9 = -0.24784313641965327321123187598392850_wp
   real(wp), parameter :: ds10 = -0.86131170741789745523421351878316690_wp
   real(wp), parameter :: ds11 = -4.2741434471979367439122664219911502E-0002_wp
   xr1 = real(h(0))
   xi1 = aimag(h(0))
   xr2 = real(h(1))+real(h(12))
   xi2 = aimag(h(1))+aimag(h(12))
   xr3 = real(h(2))+real(h(11))
   xi3 = aimag(h(2))+aimag(h(11))
   xr4 = real(h(4))+real(h(9))
   xi4 = aimag(h(4))+aimag(h(9))
   xr5 = real(h(8))+real(h(5))
   xi5 = aimag(h(8))+aimag(h(5))
   xr6 = real(h(3))+real(h(10))
   xi6 = aimag(h(3))+aimag(h(10))
   xr7 = real(h(6))+real(h(7))
   xi7 = aimag(h(6))+aimag(h(7))
   yr1 = real(h(1))-real(h(12))
   yi1 = aimag(h(1))-aimag(h(12))
   yr2 = real(h(2))-real(h(11))
   yi2 = aimag(h(2))-aimag(h(11))
   yr3 = real(h(4))-real(h(9))
   yi3 = aimag(h(4))-aimag(h(9))
   yr4 = real(h(8))-real(h(5))
   yi4 = aimag(h(8))-aimag(h(5))
   yr5 = real(h(3))-real(h(10))
   yi5 = aimag(h(3))-aimag(h(10))
   yr6 = real(h(6))-real(h(7))
   yi6 = aimag(h(6))-aimag(h(7))
   pr2 = xr2+xr5
   pi2 = xi2+xi5
   pr3 = xr3+xr6
   pi3 = xi3+xi6
   pr4 = xr4+xr7
   pi4 = xi4+xi7
   pr5 = xr2-xr5
   pi5 = xi2-xi5
   pr6 = xr3-xr6
   pi6 = xi3-xi6
   pr7 = xr4-xr7
   pi7 = xi4-xi7
   pr234 = pr2+pr3+pr4
   pi234 = pi2+pi3+pi4
   sr1 = xr1+pr234
   si1 = xi1+pi234
   ur234 = xr1+dc1*pr234
   ui234 = xi1+dc1*pi234
   ur2 = ur234+dc2*pr3+dc3*pr4
   ui2 = ui234+dc2*pi3+dc3*pi4
   ur3 = ur234+dc3*pr2+dc2*pr4
   ui3 = ui234+dc3*pi2+dc2*pi4
   ur4 = ur234+dc2*pr2+dc3*pr3
   ui4 = ui234+dc2*pi2+dc3*pi3
   ur5 = dc4*pr5+dc5*pr6+dc6*pr7
   ui5 = dc4*pi5+dc5*pi6+dc6*pi7
   ur6 = dc4*pr6+dc5*pr7-dc6*pr5
   ui6 = dc4*pi6+dc5*pi7-dc6*pi5
   ur7 = dc4*pr7-dc5*pr5-dc6*pr6
   ui7 = dc4*pi7-dc5*pi5-dc6*pi6
   sr2 = ur2+ur5
   si2 = ui2+ui5
   sr3 = ur3+ur6
   si3 = ui3+ui6
   sr4 = ur4+ur7
   si4 = ui4+ui7
   sr5 = ur2-ur5
   si5 = ui2-ui5
   sr6 = ur3-ur6
   si6 = ui3-ui6
   sr7 = ur4-ur7
   si7 = ui4-ui7

   qr1 = yr1-yr3+yr5
   qi1 = yi1-yi3+yi5
   qr3 = yr1+yr3
   qi3 = yi1+yi3
   qr_5 = yr3+yr5
   qi_5 = yi3+yi5
   qr2 = yr2-yr4+yr6
   qi2 = yi2-yi4+yi6
   qr4 = yr2+yr4
   qi4 = yi2+yi4
   qr_6 = yr4+yr6
   qi_6 = yi4+yi6
   vr1 = ds1*qr1+ds2*qr2
   vi1 = ds1*qi1+ds2*qi2
   vr2 = ds1*qr2-ds2*qr1
   vi2 = ds1*qi2-ds2*qi1
   qr34 = qr3+qr4
   qi34 = qi3+qi4
   qr_56 = qr_5+qr_6
   qi_56 = qi_5+qi_6
   vr34 = ds3*qr34-ds6*qr_56
   vi34 = ds3*qi34-ds6*qi_56
   vr3 = vr34+ds4*qr4-ds7*qr_6
   vi3 = vi34+ds4*qi4-ds7*qi_6
   vr4 = vr34+ds5*qr3-ds8*qr_5
   vi4 = vi34+ds5*qi3-ds8*qi_5
   vr56 = ds9*qr34-ds3*qr_56
   vi56 = ds9*qi34-ds3*qi_56
   vr5 = vr56+ds10*qr4-ds4*qr_6
   vi5 = vi56+ds10*qi4-ds4*qi_6
   vr6 = vr56+ds11*qr3-ds5*qr_5
   vi6 = vi56+ds11*qi3-ds5*qi_5
   tr1 = vr1+vr3
   ti1 = vi1+vi3
   tr3 = vr3-vr1-vr5
   ti3 = vi3-vi1-vi5
   tr5 = vr1-vr5
   ti5 = vi1-vi5
   tr2 = vr2+vr4
   ti2 = vi2+vi4
   tr4 = vr4-vr2-vr6
   ti4 = vi4-vi2-vi6
   tr6 = vr2-vr6
   ti6 = vi2-vi6

   h(0) = cmplx(sr1,si1,wp)
   h(1) = cmplx(sr7+ti6,si7-tr6,wp)
   h(2) = cmplx(sr6+ti5,si6-tr5,wp)
   h(3) = cmplx(sr3+ti2,si3-tr2,wp)
   h(4) = cmplx(sr5+ti4,si5-tr4,wp)
   h(5) = cmplx(sr4-ti3,si4+tr3,wp)
   h(6) = cmplx(sr2+ti1,si2-tr1,wp)
   h(7) = cmplx(sr2-ti1,si2+tr1,wp)
   h(8) = cmplx(sr4+ti3,si4-tr3,wp)
   h(9) = cmplx(sr5-ti4,si5+tr4,wp)
   h(10) = cmplx(sr3-ti2,si3+tr2,wp)
   h(11) = cmplx(sr6-ti5,si6+tr5,wp)
   h(12) = cmplx(sr7-ti6,si7+tr6,wp)
end subroutine fft13t