diff --git a/Helicity/WaveFunction/RSSpinorWaveFunction.cc b/Helicity/WaveFunction/RSSpinorWaveFunction.cc
--- a/Helicity/WaveFunction/RSSpinorWaveFunction.cc
+++ b/Helicity/WaveFunction/RSSpinorWaveFunction.cc
@@ -1,328 +1,250 @@
 // -*- C++ -*-
 //
 // RSSpinorWaveFunction.cc is a part of ThePEG - Toolkit for HEP Event Generation
 // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
 //
 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the RSSpinorWaveFunction class.
 //
 
 #include "RSSpinorWaveFunction.h"
+#include "ThePEG/Helicity/HelicityFunctions.h"
 
 using namespace ThePEG;
 using namespace ThePEG::Helicity;
 
 // calculate the Wavefunction
 void RSSpinorWaveFunction::calculateWaveFunction(unsigned int ihel) {
-  LorentzRSSpinor<double> news;
-  if(direction()==incoming) 
-    news=LorentzRSSpinor<double>(SpinorType::u);
-  else 
-    news=LorentzRSSpinor<double>(SpinorType::v);
-  unsigned int ix,iy;
-  // check helicity and type
+  // if rest frame, make sure speical case is used
+  Lorentz5Momentum ptemp=momentum();
+  double pr = ptemp.vect().mag2()/sqr(momentum().e());
+  if(pr<1e-20) {
+    ptemp.setX(ZERO);
+    ptemp.setY(ZERO);
+    ptemp.setZ(ZERO);
+  }
   assert(direction()!=intermediate);
   assert(ihel<=3);
-  // massive
   // only two valid helicities in massless case
-  assert( mass()>ZERO || ( ihel == 0 || ihel==3 ) );
-  // extract the momentum components
-  // compute the normal spinors to construct the RS spinor
-  Complex hel_wf[2][2];
-  if(direction()==incoming) {
-    // the + spinor
-    hel_wf[0][0] = 1.;
-    hel_wf[0][1] = 0.;
-    // the - spinor
-    hel_wf[1][0] = 0.;
-    hel_wf[1][1] = 1.;
-  }
-  else {
-    // the + spinor
-    hel_wf[0][0] = 0.;
-    hel_wf[0][1] = 1.;
-    // the - spinor
-    hel_wf[1][0] = 1.;
-    hel_wf[1][1] = 0.;
-  }
-  // prefactors
-  double fact = direction()==incoming ? 1. : -1.;
-  Energy pmm=mass(),pee=fact*e();
-  Energy pabs   = sqrt(sqr(px())+sqr(py())+sqr(pz()));
-  SqrtEnergy eplusp = sqrt(pee+pabs);
-  SqrtEnergy eminusp = ( pmm == ZERO ) ? ZERO : pmm/eplusp;
-  SqrtEnergy upper[2],lower[2];
-  if(direction()==incoming) {
-    upper[0] = eminusp;
-    lower[0] = eplusp ;
-    upper[1] = eplusp ;
-    lower[1] = eminusp;
-  }
-  else {
-    upper[0] =-eplusp ;
-    lower[0] = eminusp;
-    upper[1] = eminusp;
-    lower[1] =-eplusp ;
-  }
-  // now construct the spinors
-  complex<SqrtEnergy> spinor[2][4];
-  for(ix=0;ix<2;++ix) {
-    spinor[ix][0] = upper[ix]*hel_wf[ix][0];
-    spinor[ix][1] = upper[ix]*hel_wf[ix][1];
-    spinor[ix][2] = lower[ix]*hel_wf[ix][0];
-    spinor[ix][3] = lower[ix]*hel_wf[ix][1];
-  }
+  assert( mass()>ZERO || (ihel == 0 || ihel == 3 ) );
+  // new direct calculation
+  _wf = LorentzRSSpinor<double>(direction()==outgoing ? SpinorType::v : SpinorType::u);
+  // compute the spinors
+  std::array<LorentzSpinor<double>,2> spin;
+  if(ihel!=0) spin[0] = HelicityFunctions::spinor(ptemp,1,direction());
+  if(ihel!=3) spin[1] = HelicityFunctions::spinor(ptemp,0,direction());
   // compute the polarization vectors to construct the RS spinor
-  Complex vec[3][4],ii(0.,1.);
-  double ort = sqrt(0.5);
-  double r1 = ( pmm == ZERO ) ? 0. : double(pee /pmm);
-  double r2 = ( pmm == ZERO ) ? 0. : double(pabs/pmm);
-  if(direction()==incoming) {
-    vec[0][0] =-ort;
-    vec[0][1] =-ort*ii;
-    vec[0][2] = 0.;
-    vec[0][3] = 0.;
-    vec[1][0] = 0.;
-    vec[1][1] = 0.;
-    vec[1][2] = r1;
-    vec[1][3] = r2;
-    vec[2][0] = ort;
-    vec[2][1] =-ort*ii;
-    vec[2][2] = 0.;
-    vec[2][3] = 0.;
-  }
-  else {
-    vec[0][0] = ort;
-    vec[0][1] =-ort*ii;
-    vec[0][2] = 0.;
-    vec[0][3] = 0.;
-    vec[1][0] = 0.;
-    vec[1][1] = 0.;
-    vec[1][2] =-r1;
-    vec[1][3] =-r2;
-    vec[2][0] =-ort;
-    vec[2][1] =-ort*ii;
-    vec[2][2] = 0.;
-    vec[2][3] = 0.;
-  }
+  std::array<LorentzPolarizationVector,3> eps;
+  if(ihel>=2)
+    eps[0] = HelicityFunctions::polarizationVector(ptemp,2,direction(),vector_phase);
+  else
+    eps[2] = HelicityFunctions::polarizationVector(ptemp,0,direction(),vector_phase);
+  if(mass()!=ZERO&&ihel!=0&&ihel!=3)
+    eps[1] =  HelicityFunctions::polarizationVector(ptemp,1,direction());
   // now we can put the bits together to compute the RS spinor
   double or3(sqrt(1./3.)),tor3(sqrt(2./3.));
   if(ihel==3) {
-    for(ix=0;ix<4;++ix)
-      for(iy=0;iy<4;++iy)
-	news(ix,iy)=UnitRemoval::InvSqrtE*vec[0][ix]*spinor[0][iy];
+    for(unsigned int iy=0;iy<4;++iy) {
+      _wf(0,iy) = eps[0].x()*spin[0][iy];
+      _wf(1,iy) = eps[0].y()*spin[0][iy];
+      _wf(2,iy) = eps[0].z()*spin[0][iy];
+      _wf(3,iy) = eps[0].t()*spin[0][iy];
+    }
   }
   else if(ihel==2) {
-    for(ix=0;ix<4;++ix)
-      for(iy=0;iy<4;++iy)
-	news(ix,iy)=UnitRemoval::InvSqrtE*
-	  (or3*vec[0][ix]*spinor[1][iy]+tor3*vec[1][ix]*spinor[0][iy]);
+    for(unsigned int iy=0;iy<4;++iy) {
+      _wf(0,iy) = or3*eps[0].x()*spin[1][iy]+tor3*eps[1].x()*spin[0][iy];
+      _wf(1,iy) = or3*eps[0].y()*spin[1][iy]+tor3*eps[1].y()*spin[0][iy];
+      _wf(2,iy) = or3*eps[0].z()*spin[1][iy]+tor3*eps[1].z()*spin[0][iy];
+      _wf(3,iy) = or3*eps[0].t()*spin[1][iy]+tor3*eps[1].t()*spin[0][iy];
+    }
   }
   else if(ihel==1) {
-    for(ix=0;ix<4;++ix)
-      for(iy=0;iy<4;++iy)
-	news(ix,iy)=UnitRemoval::InvSqrtE*
-	  (or3*vec[2][ix]*spinor[0][iy]+tor3*vec[1][ix]*spinor[1][iy]);
+    for(unsigned int iy=0;iy<4;++iy) {
+      _wf(0,iy) = or3*eps[2].x()*spin[0][iy]+tor3*eps[1].x()*spin[1][iy];
+      _wf(1,iy) = or3*eps[2].y()*spin[0][iy]+tor3*eps[1].y()*spin[1][iy];
+      _wf(2,iy) = or3*eps[2].z()*spin[0][iy]+tor3*eps[1].z()*spin[1][iy];
+      _wf(3,iy) = or3*eps[2].t()*spin[0][iy]+tor3*eps[1].t()*spin[1][iy];
+    }
   }
   else if(ihel==0) {
-    for(ix=0;ix<4;++ix)
-      for(iy=0;iy<4;++iy)
-	news(ix,iy)=UnitRemoval::InvSqrtE*vec[2][ix]*spinor[1][iy];
-  }
-  // spinor is currently along the z axis, rotate so in right direction
-  if(pabs/pmm>1e-8) {
-    Axis axis;
-    axis.setX(fact*momentum().x()/pabs);
-    axis.setY(fact*momentum().y()/pabs);
-    axis.setZ(fact*momentum().z()/pabs);
-    double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
-    if(sinth>1e-8) {
-      LorentzRotation rot;
-      rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
-      _wf= news.transform(rot); 
+    for(unsigned int iy=0;iy<4;++iy) {
+      _wf(0,iy) = eps[2].x()*spin[1][iy];
+      _wf(1,iy) = eps[2].y()*spin[1][iy];
+      _wf(2,iy) = eps[2].z()*spin[1][iy];
+      _wf(3,iy) = eps[2].t()*spin[1][iy];
     }
-    else if (axis.z()<0.) {
-      LorentzRotation rot;
-      rot.setRotateX(Constants::pi);
-      _wf= news.transform(rot);
-    }
-    else {
-      _wf = news;
-    }    
-  }
-  else {
-    _wf=news;
   }
 }
 
 void RSSpinorWaveFunction::
 calculateWaveFunctions(vector<LorentzRSSpinor<SqrtEnergy> > & waves,
 		       tPPtr particle,Direction dir) {
   tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
     dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
   waves.resize(4);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = inspin->getProductionBasisState(ix);
     }
     else {
       inspin->decay();
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = inspin->getDecayBasisState(ix);
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<4;++ix) {
       wave.reset(ix);
       waves[ix] = wave.dimensionedWf();
     }
   }
 }
 
 void RSSpinorWaveFunction::
 calculateWaveFunctions(vector<RSSpinorWaveFunction> & waves,
 		       tPPtr particle,Direction dir) {
   tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
     dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
   waves.resize(4);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = RSSpinorWaveFunction(particle,
 					 inspin->getProductionBasisState(ix),dir);
     }
     else {
       inspin->decay();
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = RSSpinorWaveFunction(particle,
 					 inspin->getDecayBasisState(ix),dir);
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<4;++ix) {
       wave.reset(ix);
       waves[ix] = wave;
     }
   }
 }
 
 void RSSpinorWaveFunction::
 calculateWaveFunctions(vector<LorentzRSSpinor<SqrtEnergy> > & waves,
 		       RhoDMatrix & rho,
 		       tPPtr particle,Direction dir) {
   tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
     dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
   waves.resize(4);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = inspin->getProductionBasisState(ix);
       rho = RhoDMatrix(PDT::Spin3Half);
     }
     else {
       inspin->decay();
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = inspin->getDecayBasisState(ix);
       rho = inspin->rhoMatrix();
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<4;++ix) {
       wave.reset(ix);
       waves[ix] = wave.dimensionedWf();
     }
     rho = RhoDMatrix(PDT::Spin3Half);
   }
 }
 
 void RSSpinorWaveFunction::
 calculateWaveFunctions(vector<RSSpinorWaveFunction> & waves,
 		       RhoDMatrix & rho,
 		       tPPtr particle,Direction dir) {
   tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
     dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
   waves.resize(4);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = RSSpinorWaveFunction(particle,
 					 inspin->getProductionBasisState(ix),dir);
       rho = RhoDMatrix(PDT::Spin3Half);
     }
     else {
       inspin->decay();
       for(unsigned int ix=0;ix<4;++ix)
 	waves[ix] = RSSpinorWaveFunction(particle,
 					 inspin->getDecayBasisState(ix),dir);
       rho = inspin->rhoMatrix();
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<4;++ix) {
       wave.reset(ix);
       waves[ix] = wave;
     }
     rho = RhoDMatrix(PDT::Spin3Half);
   }
 }
 
 void RSSpinorWaveFunction::
 constructSpinInfo(const vector<LorentzRSSpinor<SqrtEnergy> > & waves,
 		  tPPtr particle,Direction dir,bool time) {
   assert(waves.size()==4);
   tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
     dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
   if(inspin) {
     for(unsigned int ix=0;ix<4;++ix)
       if(dir==outgoing) inspin->setBasisState(ix,waves[ix]);
       else              inspin->setDecayState(ix,waves[ix]);
   }
   else {
     RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time));
     particle->spinInfo(temp);
     for(unsigned int ix=0;ix<4;++ix)
       if(dir==outgoing) temp->setBasisState(ix,waves[ix]);
       else              temp->setDecayState(ix,waves[ix]);
   }
 }
 
 void RSSpinorWaveFunction::
 constructSpinInfo(const vector<RSSpinorWaveFunction> & waves,
 		  tPPtr particle,Direction dir,bool time) {
   assert(waves.size()==4);
   tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
     dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
   if(inspin) {
     for(unsigned int ix=0;ix<4;++ix)
       if(dir==outgoing) inspin->setBasisState(ix,waves[ix].dimensionedWf());
       else              inspin->setDecayState(ix,waves[ix].dimensionedWf());
   }
   else {
     RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time));
     particle->spinInfo(temp);
     for(unsigned int ix=0;ix<4;++ix)
       if(dir==outgoing) temp->setBasisState(ix,waves[ix].dimensionedWf());
       else              temp->setDecayState(ix,waves[ix].dimensionedWf());
   }
 }