diff --git a/Helicity/HelicityFunctions.h b/Helicity/HelicityFunctions.h
--- a/Helicity/HelicityFunctions.h
+++ b/Helicity/HelicityFunctions.h
@@ -1,172 +1,263 @@
 // -*- C++ -*-
 //
 // HelicityFunctions.h 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.
 //
 #ifndef ThePEG_HelicityFunctions_H
 #define ThePEG_HelicityFunctions_H
 //
 // This is the declaration of the HelicityFunctions header for common functions
 // used in helicity calculations to avoid duplication of code
 
 #include "ThePEG/Vectors/LorentzVector.h"
 #include "LorentzSpinor.h"
 #include "LorentzSpinorBar.h"
 
 namespace ThePEG {
 namespace Helicity {
 namespace HelicityFunctions {
 
 inline LorentzPolarizationVector polarizationVector(const Lorentz5Momentum & p,
 						    unsigned int ihel,
 						    Direction dir,
 						    VectorPhase vphase=default_vector_phase) {
   // check the direction
   assert(dir!=intermediate);
   // special helicity combination for guge invariance tests
   if(ihel==10) return p*UnitRemoval::InvE;
   // check a valid helicity combination
   assert(ihel==0 || ihel == 2 || (ihel==1 && p.mass()>ZERO));
   // convert the helicitty from 0,1,2 to -1,0,1
   int jhel=ihel-1;
   // extract the momentum components
   double fact = dir==outgoing ? -1 : 1;
   Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass();
   // calculate some kinematic quantites
   Energy2 pt2 = ppx*ppx+ppy*ppy;
   Energy pabs = sqrt(pt2+ppz*ppz);
   Energy pt = sqrt(pt2);
   // overall phase of the vector
   Complex phase(1.);
   if(vphase==vector_phase) {
     if(pt==ZERO || ihel==1) phase = 1.;
     else if(ihel==0)        phase = Complex(ppx/pt,-fact*ppy/pt);
     else                    phase = Complex(ppx/pt, fact*ppy/pt);
   }
   if(ihel!=1) phase = phase/sqrt(2.);
   // first the +/-1 helicity states
   if(ihel!=1) {
     // first the zero pt case
     if(pt==ZERO) {
       double sgnz = ppz<ZERO ? -1. : 1.;
       return LorentzPolarizationVector(-complex<double>(jhel)*phase,
 				       sgnz*phase*complex<double>(0,-fact),
 				       0.,0.);
     }
     else {
       InvEnergy opabs=1./pabs;
       InvEnergy opt  =1./pt;
       return LorentzPolarizationVector(phase*complex<double>(-jhel*ppz*ppx*opabs*opt, fact*ppy*opt),
 				       phase*complex<double>(-jhel*ppz*ppy*opabs*opt,-fact*ppx*opt),
 				       double(jhel*pt*opabs)*phase,0.);
     }
   }
   // 0 component for massive vectors
   else {
     if(pabs==ZERO) {
       return LorentzPolarizationVector(0.,0.,1.,0.);
     }
     else {
       InvEnergy empabs=pee/pmm/pabs;
       return LorentzPolarizationVector(double(empabs*ppx),double(empabs*ppy),
 				       double(empabs*ppz),double(pabs/pmm));
     }
   }
 }
 
+inline LorentzSpinor<double> spinor(const Lorentz5Momentum & p,
+				    unsigned int ihel,
+				    Direction dir) {
+  // check direction and helicity
+  assert(dir!=intermediate);
+  assert(ihel<=1);
+  // extract the momentum components
+  double fact = dir==incoming ? 1 : -1.;
+  Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass();
+  // define and calculate some kinematic quantities
+  Energy2 ptran2  = ppx*ppx+ppy*ppy;
+  Energy pabs   = sqrt(ptran2+ppz*ppz);
+  Energy ptran  = sqrt(ptran2);
+  // first need to evalulate the 2-component helicity spinors 
+  // this is the same regardless of which definition of the spinors
+  // we are using
+  complex <double> hel_wf[2];
+  // compute the + spinor for + helicty particles and - helicity antiparticles
+  if((dir==incoming && ihel== 1) || (dir==outgoing && ihel==0)) {
+    // no transverse momentum 
+    if(ptran==ZERO) {
+      if(ppz>=ZERO) {
+	hel_wf[0] = 1;
+	hel_wf[1] = 0;
+      }
+      else {
+	hel_wf[0] = 0;
+	hel_wf[1] = 1;
+      }
+    }
+    else {
+      InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
+      SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz); 
+      hel_wf[0] = denominator*rtppluspz;
+      hel_wf[1] = denominator/rtppluspz*complex<Energy>(ppx,ppy);
+    }
+  }
+  // compute the - spinor for - helicty particles and + helicity antiparticles
+  else {
+    // no transverse momentum
+    if(ptran==ZERO) {
+      if(ppz>=ZERO) {
+	hel_wf[0] = 0;
+	hel_wf[1] = 1;
+      }
+      // transverse momentum 
+      else {
+	hel_wf[0] = -1;
+	hel_wf[1] =  0;
+      }
+    }
+    else {
+      InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
+      SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
+      hel_wf[0] = denominator/rtppluspz*complex<Energy>(-ppx,ppy);
+      hel_wf[1] = denominator*rtppluspz;
+    }
+  }
+
+  SqrtEnergy upper,lower;
+  SqrtEnergy eplusp  = sqrt(max(pee+pabs,ZERO));
+  SqrtEnergy eminusp = ( pmm != ZERO ) ? pmm/eplusp : ZERO;
+  // set up the coefficients for the different cases
+  if(dir==incoming) {
+    if(ihel==1) {
+      upper = eminusp;
+      lower = eplusp;
+    }
+    else {
+      upper = eplusp;
+      lower = eminusp;
+    }
+  }
+  else {
+    if(ihel==1) {
+      upper = -eplusp;
+      lower = eminusp;
+    }
+    else {
+      upper = eminusp;
+      lower =-eplusp;
+    }
+  }
+  return LorentzSpinor<double>(upper*hel_wf[0]*UnitRemoval::InvSqrtE,
+			       upper*hel_wf[1]*UnitRemoval::InvSqrtE,
+			       lower*hel_wf[0]*UnitRemoval::InvSqrtE,
+			       lower*hel_wf[1]*UnitRemoval::InvSqrtE,
+			       (dir==incoming) ? SpinorType::u : SpinorType::v);
+}
+
+
 inline LorentzSpinorBar<double> spinorBar(const Lorentz5Momentum & p,
 					  unsigned int ihel,
 					  Direction dir) {
   // check direction and helicity
   assert(dir!=intermediate);
   assert(ihel<=1);
   // extract the momentum components
   double fact = dir==incoming ? 1. : -1.;
   Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass();
   // define and calculate some kinematic quantities
   Energy2 ptran2  = ppx*ppx+ppy*ppy;
   Energy pabs   = sqrt(ptran2+ppz*ppz);
   Energy ptran  = sqrt(ptran2);
   // first need to evalulate the 2-component helicity spinors
   Complex hel_wf[2];
   // compute the + spinor for + helicty particles and - helicity antiparticles
   if((dir==outgoing && ihel== 1) || (dir==incoming && ihel==0)) {
     // no transverse momentum
     if(ptran==ZERO) {
       if(ppz>=ZERO) {
   	hel_wf[0] = 1;
   	hel_wf[1] = 0;
       }
       else {
   	hel_wf[0] = 0;
   	hel_wf[1] = 1;
       }
     }
     else {
       InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
       SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
       hel_wf[0] = denominator*rtppluspz;
       hel_wf[1] = denominator/rtppluspz*complex<Energy>(ppx,-ppy);
     }
   }
   // compute the - spinor for - helicty particles and + helicity antiparticles
   else {
     // no transverse momentum
     if(ptran==ZERO) {
       if(ppz>=ZERO) {
   	hel_wf[0] = 0;
   	hel_wf[1] = 1;
       }
       // transverse momentum
       else {
   	hel_wf[0] = -1;
   	hel_wf[1] =  0;
       }
     }
     else {
       InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
       SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
       hel_wf[0] = denominator/rtppluspz*complex<Energy>(-ppx,-ppy);
       hel_wf[1] = denominator*rtppluspz;
     }
   }
   SqrtEnergy upper, lower;
   SqrtEnergy eplusp  = sqrt(max(pee+pabs,ZERO));
   SqrtEnergy eminusp = ( pmm!=ZERO ) ? pmm/eplusp : ZERO;
   // set up the coefficients for the different cases
   if(dir==outgoing) {
     if(ihel==1) {
       upper = eplusp;
       lower = eminusp;
     }
     else {
       upper = eminusp;
       lower = eplusp;
     }
   }
   else {
     if(ihel==1) {
   	upper = eminusp;
   	lower = -eplusp;
     }
     else {
       upper =-eplusp;
       lower = eminusp;
     }
   }
   // now finally we can construct the spinors
   return LorentzSpinorBar<double>(upper*hel_wf[0]*UnitRemoval::InvSqrtE,
 				  upper*hel_wf[1]*UnitRemoval::InvSqrtE,
 				  lower*hel_wf[0]*UnitRemoval::InvSqrtE,
 				  lower*hel_wf[1]*UnitRemoval::InvSqrtE,
 				  (dir==incoming) ? SpinorType::v : SpinorType::u);
 }
 
 }
 }
 }
   
 #endif /* ThePEG_HelicityFunctions_H */
diff --git a/Helicity/WaveFunction/SpinorWaveFunction.cc b/Helicity/WaveFunction/SpinorWaveFunction.cc
--- a/Helicity/WaveFunction/SpinorWaveFunction.cc
+++ b/Helicity/WaveFunction/SpinorWaveFunction.cc
@@ -1,340 +1,250 @@
 // -*- C++ -*-
 //
 // SpinorWaveFunction.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 SpinorWaveFunction class.
 //
 // Author: Peter Richardson
 //
 
 #include "SpinorWaveFunction.h"
 #include "ThePEG/Repository/CurrentGenerator.h"
 #include "SpinorBarWaveFunction.h"
+#include "ThePEG/Helicity/HelicityFunctions.h"
 
 using namespace ThePEG;
 using namespace Helicity;
 
 // calculate the Wavefunction
 void SpinorWaveFunction::calculateWaveFunction(unsigned int ihel) {
-  // check helicity is O.K.
-  Direction dir = direction();
-  if(dir==intermediate) throw ThePEG::Helicity::HelicityConsistencyError() 
-    << "In SpinorWaveFunction::calcluateWaveFunction "
-    << "particle must be incoming or outgoing not intermediate" 
-    << Exception::abortnow;
-  if(ihel>1) throw ThePEG::Helicity::HelicityConsistencyError() 
-    << "Invalid Helicity = " << ihel << " requested for Spinor" 
-    << Exception::abortnow;
-  // extract the momentum components
-  double fact=-1.; if(dir==incoming){fact=1.;}
-  Energy ppx=fact*px(),ppy=fact*py(),ppz=fact*pz(),pee=fact*e(),pmm=mass();
-  // define and calculate some kinematic quantities
-  Energy2 ptran2  = ppx*ppx+ppy*ppy;
-  Energy pabs   = sqrt(ptran2+ppz*ppz);
-  Energy ptran  = sqrt(ptran2);
-  // first need to evalulate the 2-component helicity spinors 
-  // this is the same regardless of which definition of the spinors
-  // we are using
-  complex <double> hel_wf[2];
-  // compute the + spinor for + helicty particles and - helicity antiparticles
-  if((dir==incoming && ihel== 1) || (dir==outgoing && ihel==0)) {
-    // no transverse momentum 
-    if(ptran==ZERO) {
-      if(ppz>=ZERO) {
-	hel_wf[0] = 1;
-	hel_wf[1] = 0;
-      }
-      else {
-	hel_wf[0] = 0;
-	hel_wf[1] = 1;
-      }
-    }
-    else {
-      InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
-      SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz); 
-      hel_wf[0] = denominator*rtppluspz;
-      hel_wf[1] = denominator/rtppluspz*complex<Energy>(ppx,ppy);
-    }
-  }
-  // compute the - spinor for - helicty particles and + helicity antiparticles
-  else {
-    // no transverse momentum
-    if(ptran==ZERO) {
-      if(ppz>=ZERO) {
-	hel_wf[0] = 0;
-	hel_wf[1] = 1;
-      }
-      // transverse momentum 
-      else {
-	hel_wf[0] = -1;
-	hel_wf[1] =  0;
-      }
-    }
-    else {
-      InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
-      SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
-      hel_wf[0] = denominator/rtppluspz*complex<Energy>(-ppx,ppy);
-      hel_wf[1] = denominator*rtppluspz;
-    }
-  }
-
-  SqrtEnergy upper,lower;
-  SqrtEnergy eplusp  = sqrt(max(pee+pabs,ZERO));
-  SqrtEnergy eminusp = ( pmm != ZERO ) ? pmm/eplusp : ZERO;
-  // set up the coefficients for the different cases
-  if(dir==incoming) {
-    if(ihel==1) {
-      upper = eminusp;
-      lower = eplusp;
-    }
-    else {
-      upper = eplusp;
-      lower = eminusp;
-    }
-  }
-  else {
-    if(ihel==1) {
-      upper = -eplusp;
-      lower = eminusp;
-    }
-    else {
-      upper = eminusp;
-      lower =-eplusp;
-    }
-  }
-  // now finally we can construct the spinors
-  _wf = LorentzSpinor<double>( (dir==incoming) ? SpinorType::u : SpinorType::v);
-  _wf[0]=upper*hel_wf[0]*UnitRemoval::InvSqrtE;
-  _wf[1]=upper*hel_wf[1]*UnitRemoval::InvSqrtE;
-  _wf[2]=lower*hel_wf[0]*UnitRemoval::InvSqrtE;
-  _wf[3]=lower*hel_wf[1]*UnitRemoval::InvSqrtE;
+  _wf = HelicityFunctions::spinor(momentum(),ihel,direction());
 }
 
 SpinorBarWaveFunction SpinorWaveFunction::bar() {
   Lorentz5Momentum p = momentum();
   if(direction()==outgoing) p *= -1.;
   tcPDPtr ptemp = particle();
   if(direction()==incoming&&particle()->CC())
     ptemp = particle()->CC();
   return SpinorBarWaveFunction(p,ptemp,_wf.bar(),direction());
 }	  
 
 void SpinorWaveFunction::
 calculateWaveFunctions(vector<LorentzSpinor<SqrtEnergy> > & waves,
 		       tPPtr particle,Direction dir) {
   tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
     dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
   waves.resize(2);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       waves[0] = inspin->getProductionBasisState(0);
       waves[1] = inspin->getProductionBasisState(1);
     }
     else {
       inspin->decay();
       if( (particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) || 
 	  (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
 	waves[0] = inspin->getDecayBasisState(0).conjugate();
 	waves[1] = inspin->getDecayBasisState(1).conjugate();
       }
       else {
 	waves[0] = inspin->getDecayBasisState(0);
 	waves[1] = inspin->getDecayBasisState(1);
       }
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<2;++ix) {
       wave.reset(ix);
       waves[ix] = wave.dimensionedWave();
     }
   }
 }
 
 void SpinorWaveFunction::
 calculateWaveFunctions(vector<SpinorWaveFunction> & waves,
 		       tPPtr particle,Direction dir) {
   tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
     dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
   waves.resize(2);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       for(unsigned int ix=0;ix<2;++ix)
 	waves[ix] = SpinorWaveFunction(particle,
 				       inspin->getProductionBasisState(ix),dir);
     }
     else {
       inspin->decay();
       if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
 	 (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
 	for(unsigned int ix=0;ix<2;++ix)
 	  waves[ix] = SpinorWaveFunction(particle,
 					 inspin->getDecayBasisState(ix).conjugate(),dir);
       }
       else {
 	for(unsigned int ix=0;ix<2;++ix)
 	  waves[ix] = SpinorWaveFunction(particle,
 					 inspin->getDecayBasisState(ix),dir);
       }
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<2;++ix) {
       wave.reset(ix);
       waves[ix] = wave;
     }
   }
 }
 
 void SpinorWaveFunction::
 calculateWaveFunctions(vector<LorentzSpinor<SqrtEnergy> > & waves,
 		       RhoDMatrix & rho,
 		       tPPtr particle,Direction dir) {
   tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
     dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
   waves.resize(2);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       waves[0] = inspin->getProductionBasisState(0);
       waves[1] = inspin->getProductionBasisState(1);
       rho = RhoDMatrix(PDT::Spin1Half);
     }
     else {
       inspin->decay();
       if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
 	 (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
 	waves[0] = inspin->getDecayBasisState(0).conjugate();
 	waves[1] = inspin->getDecayBasisState(1).conjugate();
       }
       else {
 	waves[0] = inspin->getDecayBasisState(0);
 	waves[1] = inspin->getDecayBasisState(1);
       }
       rho = inspin->rhoMatrix();
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<2;++ix) {
       wave.reset(ix);
       waves[ix] = wave.dimensionedWave();
     }
     rho = RhoDMatrix(PDT::Spin1Half);
   }
 }
 
 void SpinorWaveFunction::
 calculateWaveFunctions(vector<SpinorWaveFunction> & waves,
 		       RhoDMatrix & rho,
 		       tPPtr particle,Direction dir) {
   tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
     dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
   waves.resize(2);
   // spin info object exists
   if(inspin) {
     if(dir==outgoing) {
       for(unsigned int ix=0;ix<2;++ix)
 	waves[ix] = SpinorWaveFunction(particle,
 				       inspin->getProductionBasisState(ix),dir);
       rho = RhoDMatrix(PDT::Spin1Half);
     }
     else {
       inspin->decay();
       if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
 	 (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
 	for(unsigned int ix=0;ix<2;++ix)
 	  waves[ix] = SpinorWaveFunction(particle,
 					 inspin->getDecayBasisState(ix).conjugate(),dir);
       }
       else {
 	for(unsigned int ix=0;ix<2;++ix)
 	  waves[ix] = SpinorWaveFunction(particle,
 					 inspin->getDecayBasisState(ix),dir);
       }
       rho = inspin->rhoMatrix();
     }
   }
   // do the calculation
   else {
     assert(!particle->spinInfo());
     SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
     for(unsigned int ix=0;ix<2;++ix) {
       wave.reset(ix);
       waves[ix] = wave;
     }
     rho = RhoDMatrix(PDT::Spin1Half);
   }
 }
 
 void SpinorWaveFunction::
 constructSpinInfo(const vector<LorentzSpinor<SqrtEnergy> > & waves,
 		  tPPtr particle,Direction dir,bool time) {
   assert(waves.size()==2);
   tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
     dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
   if(inspin) {
     for(unsigned int ix=0;ix<2;++ix) {
       if(( dir == outgoing &&  time) || 
 	 ( dir == incoming && !time))
 	inspin->setBasisState(ix,waves[ix]);
       else
 	inspin->setDecayState(ix,waves[ix]);
     }
   }
   else {
     FermionSpinPtr temp = new_ptr(FermionSpinInfo(particle->momentum(),time));
     particle->spinInfo(temp);
     for(unsigned int ix=0;ix<2;++ix) {
       if(( dir == outgoing &&  time) || 
 	 ( dir == incoming && !time))
 	temp->setBasisState(ix,waves[ix]);
       else
 	temp->setDecayState(ix,waves[ix]);
     }
   }
 }
 
 void SpinorWaveFunction::
 constructSpinInfo(const vector<SpinorWaveFunction> & waves,
 		  tPPtr particle,Direction dir,bool time) {
   assert(waves.size()==2);
   tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
     dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
   if(inspin) {
     for(unsigned int ix=0;ix<2;++ix) {
       if(( dir == outgoing &&  time) || 
 	 ( dir == incoming && !time)) 
 	inspin->setBasisState(ix,waves[ix].dimensionedWf());
       else
 	inspin->setDecayState(ix,waves[ix].dimensionedWf());
     }
   }
   else {
     FermionSpinPtr temp = new_ptr(FermionSpinInfo(particle->momentum(),time));
     particle->spinInfo(temp);
     for(unsigned int ix=0;ix<2;++ix) {
       if(( dir == outgoing &&  time) || 
 	 ( dir == incoming && !time))
 	temp->setBasisState(ix,waves[ix].dimensionedWf());
       else
 	temp->setDecayState(ix,waves[ix].dimensionedWf());
     }
   }
 }