static unsigned int is_num_laser(particle* particles) {
    singleton_geometry *geometry = get_singleton_geometry();
    unsigned int n = geometry->n;

    double x_foc = 50.e-6;
    double y_foc = 100.e-6;
    double rb = 25.e-06;
    unsigned int iter = 0;

    for (unsigned int i = 0; i < n; i++) {
        if(particles[i].bnd || particles[i].molten || particles[i].states==INACTIVE) continue;
        double xp = particles[i].px;
        double yp = particles[i].py;
        double r2 = (xp-x_foc)*(xp-x_foc) + (yp-y_foc)*(yp-y_foc);
        double dist = sqrt(r2);

        /*
           99% of beam's power will flow through a circle of radius r = 1.52 ⋅ w ( z )
           {\displaystyle r=1.52\cdot w(z)} r=1.52\cdot w(z).[wiki: Gaussian Beam]
        */

        if(particles[i].states==LASER && dist<=1.52*rb) iter++;
        //if(particles[i].states==LASER && dist<=rb) iter++;
        //if(particles[i].states==LASER) iter++;
    }
    return iter;
}