subroutine Box_Muller_polar(y1,y2)

       real :: y1,y2
       real :: x1,x2,w
       real,dimension(2) :: u
      
       w=2.0

       do while( w >= 1.0 ) ;
       call random_number(u(1:2))
           x1 = 2.0 * u(1) - 1.0;
           x2 = 2.0 * u(2) - 1.0;
           w = x1 * x1 + x2 * x2;
       end do

         w = sqrt( (-2.0 * log( w ) ) / w );
         y1 = x1 * w;
         y2 = x2 * w;

       end subroutine Box_Muller_polar