anchor.diff:  
     if(cdist_ampX < pieni) then
       do j=1,napx
  -      phiX   = twopi*real(rndm4(),fPrec)
  -      iiX    = (-one*cdist_emitX) * log_mb(real(rndm4(),fPrec))
  +      phiX   = twopi*rndm4()
  +      iiX    = (-one*cdist_emitX) * log_mb(rndm4())
         xv1(j) = sqrt((two*iiX)*cdist_betaX) * cos_mb(phiX)
         yv1(j) = (-one*sqrt((two*iiX)/cdist_betaX)) * (sin_mb(phiX) + cdist_alphaX * cos_mb(phiX))
       end do
  
positive.code: 
    if(cdist_ampY > zero) then
      do j=1, napx
        emitY  = cdist_emitY*(cdist_ampY + ((two*real(rndm4()-half,fPrec))*cdist_smearY))**2
        sigmaY = sqrt(cdist_betaY*emitY)
        xv2(j) = sigmaY * sin_mb(twopi*real(rndm4(),fPrec))
        if(rndm4() > half) then
          yv2(j) = sqrt(emitY/cdist_betaY-xv2(j)**2/cdist_betaY**2)-(cdist_alphaY*xv2(j))/cdist_betaY
        else
  
negative.code: 
    if(cdist_ampY < pieni) then
      do j=1,napx
        phiY   = twopi*rndm4()
        iiY    = (-one*cdist_emitY) * log_mb(rndm4())
        xv2(j) = sqrt((two*iiY)*cdist_betaY) * cos_mb(phiY)
        yv2(j) = (-one*sqrt((two*iiY)/cdist_betaY)) * (sin_mb(phiY) + cdist_alphaY * cos_mb(phiY))
      end do
  