Note
This subroutine initializes a real matrix tab of dimensions nx by ny with ones
subroutine make_scratches() !================================================================================================ !< @note This subroutine initializes a real matrix `tab` of dimensions `nx` by `ny` with ones ! and adds `sn` linear scratches with value 0. Each scratch has a random length up to `sl`, ! a constant width `sw`, and a random orientation. Scratches can touch the matrix boundaries. ! ! @endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4) :: nx ! Number of rows integer(kind=I4) :: ny ! Number of columns integer(kind=I4) :: sn ! Number of scratches integer(kind=I4) :: sw ! Scratch width integer(kind=I4) :: sl ! Maximum scratch length integer(kind=I4) :: i, j, n real(kind=R8) :: theta, length, depth, x0, y0, x1, y1 real(kind=R8) :: rand, r1, r2 nx = PARAM%width ny = PARAM%height read(JOB,*) sn, sw, sl ; LINE_READ = LINE_READ + 1 ; write(SPY,*) "line: ", LINE_READ, "sn, sw, sl ", sn, sw, sl ! Initialize matrix with ones PARAM%surf(1:nx, 1:ny) = 1.0_R8 ! Seed random number generator call random_seed() ! Generate each scratch do n = 1, sn ! Random starting point call random_number(rand) x0 = rand * nx call random_number(rand) y0 = rand * ny ! Random depth ! Génération de deux variables gamma call gamma_random(shape = 3._R8, output = r1) call gamma_random(shape = 1._R8, output = r2) ! Transformation pour obtenir la loi bêta depth = r1 / (r1 + r2) ! Random orientation (angle in radians) call random_number(rand) theta = rand * 2.0 * PI_R8 ! Random length up to sl call random_number(rand) length = ( 2 * (rand - 0.5_R8) / 10 + 0.9_R8 ) * sl ! Calculate end point of the scratch x1 = x0 + length * cos(theta) y1 = y0 + length * sin(theta) ! Draw the scratch with width sw do i = max(1, floor(min(x0, x1)) - sw), min(nx, ceiling(max(x0, x1)) + sw) do j = max(1, floor(min(y0, y1)) - sw), min(ny, ceiling(max(y0, y1)) + sw) ! Check if point (i,j) is within sw/2 distance from the line segment if ( point_to_line_distance(real(i, kind=R8), real(j, kind=R8), x0, y0, x1, y1) <= sw/2.0 ) then PARAM%surf(i, j) = depth endif enddo enddo enddo contains subroutine gamma_random(shape, output) implicit none real(kind=R8), intent( in) :: shape real(kind=R8), intent(out) :: output real(kind=R8) :: b, c real(kind=R8) :: u, v, w, y ! algorithme de marsaglia pour générer une variable gamma b = shape - 1.0/3.0 c = 1.0 / sqrt(9.0 * b) do call random_number(u) call random_number(v) v = 1.0 + c * log(v) / sqrt(u) if (v <= 0.0) cycle w = v * v * v call random_number(y) if (y < 1.0 - 0.0331 * (log(v) ** 4)) exit if (log(y) < 0.5 * log(v) * log(v) + b * (1.0 - v + log(v))) exit enddo output = b * w return endsubroutine gamma_random !----------------------------------------- real(kind=R8) function point_to_line_distance(px, py, x0, y0, x1, y1) !! Calculate the shortest distance from a point to a line segment. implicit none real(kind=R8), intent(in) :: px !! x-coordinate of the point real(kind=R8), intent(in) :: py !! y-coordinate of the point real(kind=R8), intent(in) :: x0 !! x-coordinate of the line segment start real(kind=R8), intent(in) :: y0 !! y-coordinate of the line segment start real(kind=R8), intent(in) :: x1 !! x-coordinate of the line segment end real(kind=R8), intent(in) :: y1 !! y-coordinate of the line segment end real(kind=R8) :: dx, dy, t, proj_x, proj_y dx = x1 - x0 dy = y1 - y0 ! Handle case where start and end points are the same if ( (dx**2 + dy**2) < 1.e-12 ) then point_to_line_distance = sqrt((px - x0)**2 + (py - y0)**2) return endif ! Project point onto the line t = max(0.0_R8, min(1.0_R8, ((px - x0)*dx + (py - y0)*dy)/(dx*dx + dy*dy))) proj_x = x0 + t * dx proj_y = y0 + t * dy ! Calculate distance point_to_line_distance = sqrt( (px - proj_x)**2 + (py - proj_y)**2 ) return endfunction point_to_line_distance !----------------------------------------- endsubroutine make_scratches