    {
       posX[threadIdx.x] = (float)(DX*injectWidth*randPool[randOffset+6]+botXStart);
       posY[threadIdx.x] = (float)(DY*randPool[randOffset+19]);
-      vpar = (float)((1.414f*rsqrtf(D_SIGMA1*velmass))*
+      vpar = (float)((1.414f*rsqrtf(D_SIGMA_CI*velmass))*
          sqrtf(-logf(1.0f-randPool[randOffset+20] + FLT_MIN)));
       tpar = (float)(D_TPI*randPool[randOffset+21] - D_PI);
       velX[threadIdx.x] = (float)vpar*__sinf((float)tpar);
-      vpar = (float)((1.414f*rsqrtf(D_SIGMA1*velmass))*
+      vpar = (float)((1.414f*rsqrtf(D_SIGMA_CI*velmass))*
          sqrtf(-logf(1.0f-randPool[randOffset+22] + FLT_MIN)));
+      // For sincos I need a range of -pi to pi
       tpar = (float)(D_TPI*randPool[randOffset+23] - D_PI);
       __sincosf(tpar, &stpar, &ctpar);
-      velY[threadIdx.x] = vpar*stpar;
-      velZ[threadIdx.x] = vpar*ctpar;
-      posY[threadIdx.x] = posY[threadIdx.x]+D_DELT*velY[threadIdx.x];
-
+      velY[threadIdx.x] = vpar*stpar + (1.1f*rsqrtf(D_SIGMA_CI*velmass));
+      velZ[threadIdx.x] = vpar*ctpar;
+      posY[threadIdx.x] = posY[threadIdx.x]+D_DELT*velY[threadIdx.x];
       posY[threadIdx.x] = max(posY[threadIdx.x], 0.0f);
+
       injectWriteBlock(ionColdLoc, ionColdVel, numIonsCold, 
          posX, posY, velX, velY, velZ);
    }

   {
      posX[threadIdx.x] = (float)(DX*injectWidth*randPool[randOffset+6]+botXStart);
      posY[threadIdx.x] = (float)(DY*randPool[randOffset+7]);
      vpar = (float)((1.414f*rsqrtf(D_SIGMA))*
         sqrtf(-logf(1-randPool[randOffset+8] + FLT_MIN)));
      tpar = (float)(D_TPI*randPool[randOffset+9] - D_PI);
      velX[threadIdx.x] = (float)(vpar*__sinf(tpar));
      vpar = (float)((1.414f/sqrtf(D_SIGMA))*
         sqrtf(-logf(1-randPool[randOffset+10] + FLT_MIN)));
      tpar = (float)(D_TPI*randPool[randOffset+11] - D_PI);
      __sincosf(tpar, &stpar, &ctpar);
      velY[threadIdx.x] = (float)(vpar*stpar);
      velY[threadIdx.x] = static_cast<float>(velY[threadIdx.x]+(1.1f*rsqrtf(D_SIGMA)));
      velZ[threadIdx.x] = vpar*ctpar;
      posY[threadIdx.x] = posY[threadIdx.x]+D_DELT*velY[threadIdx.x];

      posY[threadIdx.x] = max(posY[threadIdx.x], 0.0f);
      injectWriteBlock(eleColdLoc, eleColdVel, numElectronsCold, 
         posX, posY, velX, velY, velZ);
   }

   {
      posX[threadIdx.x] = (float)(DX*NX1*randPool[randOffset]);
      posY[threadIdx.x] = (float)(DY*(NY1-1)+DY*randPool[randOffset+1]);
      vpar=(float)((1.414f*rsqrtf(D_SIGMA_HE))*
         sqrtf(-logf(1.0f-randPool[randOffset+2] + FLT_MIN)));
      tpar = (float)(D_TPI*randPool[randOffset+3] - D_PI);
      velX[threadIdx.x] = (float)vpar*__sinf((float)tpar);
      vpar=(float)((1.414f*rsqrtf(D_SIGMA_HE))*
         sqrtf(-logf(1.0f-randPool[randOffset+4] + FLT_MIN)));
      // For sincos I need a range of -pi to pi
      tpar=(float)(D_TPI*randPool[randOffset+5] - D_PI);
      __sincosf(tpar, &stpar, &ctpar);
      velY[threadIdx.x] = vpar*stpar - (1.1f*rsqrtf(D_SIGMA_HE));
      velZ[threadIdx.x] = vpar*ctpar;
      posY[threadIdx.x] = posY[threadIdx.x]+D_DELT*velY[threadIdx.x];

