while (z > -2.5) // z: posiotion, 2.5 = hight of fall
{
    t = t + dt; // t: time, dt: time step

    double sfz = 0, staw = 0; //sfz: total vertical force, stw: total vertical torque
    
    for(double l = 0 ; l < L ; l = l + dl) //integrating the force from l = 0 to l = L, L: length of blades
    {
        double dfz, dtaw, alpha, beta, dA, dfl, dfd, Cl, Cd; //all these parameters will be defined
				
        beta = atan(-vz / (om * l)); //beta: angle of air flow with horizon 
        alpha = beta - tet; //alpha: angle of attack, tet: angle of blades with horizon
                
        dA = w * dl; //dA = area of an infitesimal segment along the blade, w = width of the blade
                
        Cl =  l_coe[(int)(2 * alpha * 180 / 3.14159265)]; //reading the lift and drag coefficients as a function of 
        Cd =  d_coe[(int)(2 * alpha * 180 / 3.14159265)]; //angle of attack (alpha) from an external file resulted from FLUENT
		
        dfd = rho / 2 * dA * Cd * (vz*vz + l*l*om*om); //infinitesimal drag force, applied to a small portion of the blade 
        dfl = rho / 2 * dA * Cl * (vz*vz + l*l*om*om); //infinitesimal lift force, applied to a small portion of the blade 
		
        dfz = dfl * cos(beta) + dfd * sin(beta); //vertical components of drag and lift forces
        dtaw = (dfl * sin(beta) - dfd * cos(beta))*l; //vertical torque resulted from drag and lift forces
		
        sfz = sfz + dfz; //summing up the infinitesimal vertical forces sum(n+1) = sum(n) + df
        staw = staw + dtaw; //summing up the infinitesimal vertical torques
    }
            
    
    sfz = sfz * n; //n: number of blades
    staw = staw * n;


    
    az = sfz / m - g; //az: linear acceleration, g: gravity accelration, m: mass of the device
    alp = staw / I; //alp: angular acceleration, I: rotational inertia
            
    vz = vz + az * dt; //vz: vertical velocity
    om = om + alp * dt; //om: angular velocity
        
    z = z + vz * dt + 0.5 * az * dt * dt; //z: position
}