Note
Function that calculates the mean absolute difference between the desired Acf and the one obtained. However, the important zone where both should match is above acf__z - where the correlation lengths are determined.
If the mean absolute difference is below the criterion, the loops to improve the acf are stopped.
The function can also plot the acfs.
subroutine plt__acf() !================================================================================================ !<@note Function that calculates the mean absolute difference between the desired Acf and !< the one obtained. !< However, the important zone where both should match is above acf__z - where the correlation !< lengths are determined. !< !< If the mean absolute difference is below the criterion, the loops to improve the acf are !< stopped. !< !< The function can also plot the acfs. !< !<@endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4) :: i, w, h, mw, mh, uplot, ll logical(kind=I4) :: is_x, is_y real (kind=R8) :: lim_crit_acf, crit_acf, dxy, l1, l2 character(len=512) :: plt_acf real(kind=R8), allocatable, dimension(:,:) :: tab_tmp real(kind=R8), dimension(1:8) :: res read(JOB,*) lim_crit_acf ; LINE_READ = LINE_READ +1 ; write(SPY,*) LINE_READ, "Acf criterion ", lim_crit_acf read(JOB,*) plt_acf ; LINE_READ = LINE_READ +1 ; write(SPY,*) LINE_READ, trim(plt_acf) PARAM%crt_acf = lim_crit_acf ! mean absolute difference limit w = PARAM%width h = PARAM%height ! Mean absolute difference between the calculated acf and the theoretical acf !................................................................................... allocate( tab_tmp(1:w, 1:h) ) tab_tmp(1:w, 1:h) = 0 where( PARAM%imp_acf > PARAM%acf__z ) tab_tmp = abs( PARAM%acf_surf - PARAM%imp_acf ) crit_acf = 100 * sum( tab_tmp(1:w, 1:h) ) / count( tab_tmp(1:w, 1:h) > 0 ) PARAM%res_acf = crit_acf deallocate( tab_tmp ) write(TER,*) "acf difference ", crit_acf write(SPY,*) 'acf difference ', crit_acf ! if the acf criterion is reached, the loops stop: a means is to modify the max number ! of loops, so that the main loop is exited. if (lim_crit_acf > 0. .and. lim_crit_acf > crit_acf) NB_ITER = 1 !................................................................................... ! Graphs? !................................................................................... ! if 'x' is present, plot the graph along the principal axis ! if 'y' is present, plot the graph along the secondary axis is_x = ( index(trim(plt_acf), 'x') /= 0 ) is_y = ( index(trim(plt_acf), 'y') /= 0 ) if ( .not.( is_x .or. is_y ) ) return ! ensure that w and h are odd, or act accordingly mw = w / 2 mh = h / 2 if ( w == 2 * (w/2) ) mw = w/2 + 1 if ( h == 2 * (h/2) ) mh = h/2 + 1 call ellipse_acf( tabin = PARAM%acf_surf(1:w, 1:h), & ! IN long = w, & ! IN larg = h, & ! IN p_acv = res(1:8), & ! OUT -> correlation lengths cut = PARAM%acf__z, & ! IN -> z cut plane scale_xy = [PARAM%surf_dx, PARAM%surf_dy], & ! IN -> lags along x and y omp = .true. ) ! IN -> use multithread? write(TER,*) res(1:2), res(4) write(SPY,*) 'acf lengths and roughness orientation ', res(1:2), res(4) ! parameters for the plot ll = 2 * min(mw, mh) - 3 dxy = sqrt( PARAM%surf_dx**2 + PARAM%surf_dy**2 ) call get_unit( uplot ) if ( is_x ) call graph(axis = 1) if ( is_y ) call graph(axis = 2) contains subroutine graph(axis) !================================================================================================ !<@note Function that plots the graphs to compare the ACF along the primary and/or secondary axes. !<@endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4), intent(in) :: axis !! *1 or 2 for primary or secondary axis* character(len=256) :: file_acf, file_gpl, title character(len=003) :: xlab real(kind=R8) :: angle real(kind=R8), allocatable, dimension(:) :: profile_acf_surf, profile_imp_acf allocate( profile_acf_surf(1:ll) ) allocate( profile_imp_acf (1:ll) ) if ( axis == 1 ) then angle = PARAM%a_acf file_acf = 'out/gpl/acfx.txt' file_gpl = 'out/gpl/acfx.gpl' title = '"ACF comparison along primary axis X"' xlab = '"X"' else angle = PARAM%a_acf + 90. file_acf = 'out/gpl/acfy.txt' file_gpl = 'out/gpl/acfy.gpl' title = '"ACF comparison along secondary axis Y"' xlab = '"Y"' endif ! extract the ACF profile along a particular direction call profile_at_angle( tab = PARAM%acf_surf(1:w, 1:h), profile = profile_acf_surf(1:ll), theta = angle ) call profile_at_angle( tab = PARAM%imp_acf (1:w, 1:h), profile = profile_imp_acf (1:ll), theta = angle ) open(uplot, file = trim(file_acf)) write(uplot, *) 'X', '"calculated acf"', '"theoretical acf"' do i = 1, ll write(uplot, *) (i - ll/2) * dxy, real(profile_acf_surf(i), kind=R4), & ! real(profile_imp_acf (i), kind=R4) ! if ( i - ll/2 < 0 ) then if ( profile_acf_surf(i) < PARAM%acf__z .and. profile_acf_surf(i + 1) > PARAM%acf__z ) l1 = (i - ll/2) * dxy endif if ( i - ll/2 > 0 .and. i < ll ) then if ( profile_acf_surf(i) > PARAM%acf__z .and. profile_acf_surf(i + 1) < PARAM%acf__z ) l2 = (i - ll/2) * dxy endif enddo close(uplot) open(uplot, file = trim(file_gpl)) write(uplot, '(a)') 'set title ' // trim(title) write(uplot, '(a)') 'set xlabel ' // trim(xlab) write(uplot, '(a)') 'set ylabel "ACF"' write(uplot, '(a,f4.2,a,f5.2,a)') "set arrow from graph 0, first ", PARAM%acf__z, & ! " to graph 1, first ", PARAM%acf__z, ' nohead lc rgb "#000000" front' ! write(uplot, '(a,E8.2,a,E8.2,a,f5.2,a)') "set arrow from ", l1, ", graph 0 to ", & ! l1, ",", PARAM%acf__z, ' nohead lc rgb "#000000" front' ! write(uplot, '(a,E8.2,a,E8.2,a,f5.2,a)') "set arrow from ", l2, ", graph 0 to ", & ! l2, ",", PARAM%acf__z, ' nohead lc rgb "#000000" front' ! write(uplot, '(a,E8.2,a,E8.2,a,f5.2)') 'set label "L1 = ', res(axis), '" at ', l2, ',', PARAM%acf__z + 0.1 write(uplot, '(a,i2,a)' ) 'plot "' // trim(file_acf) // '" using 1:2 with lines title "acf real surface", "' // trim(file_acf) // '" using 1:3 with lines title "theoretical acf"' write(uplot, '(a)' ) 'pause -1 "Hit return to continue"' write(uplot, '(a)' ) 'q' close(uplot) call system('gnuplot "' // trim(file_gpl) // '"') deallocate( profile_acf_surf ) deallocate( profile_imp_acf ) return endsubroutine graph subroutine profile_at_angle(tab, profile, theta) !================================================================================================ !<@note Function that extract the ACF profile along a particular direction !<@endnote !------------------------------------------------------------------------------------------------ implicit none real(kind=R8), intent(in), dimension(1:w, 1:h) :: tab real(kind=R8), intent(in) :: theta real(kind=R8), intent(out), dimension(1:ll) :: profile integer(kind=I4) :: p, nx, ny real (kind=R8) :: r, x, y, xb, yb, xm, ym, xp, yp, h1, h2, h3, h4, hh do p = -ll / 2, ll / 2 ! identifying a point on the diameter r = p ! corresponding algebraic radius ! projection on x and y of the point marked by its radius and angle ! by taking the lower integer, we have the number of the bottom row and left-hand column of the rectangle ! the remainder (x-nx) represents the abscissa of the point in the rectangle with sides=1 ! the 0.9999 coefficient is used to avoid falling right on an existing point x = mw + r * cos(theta * PI_R8 / 180) * 0.9999_R8 ; nx = floor(x) ; xb = x -nx y = mh + r * sin(theta * PI_R8 / 180) * 0.9999_R8 ; ny = floor(y) ; yb = y -ny xm = UN -xb ; xp = xb ym = UN -yb ; yp = yb if ( nx+1 <= w .and. ny+1 <= h .and. & ! nx >= 1 .and. ny >= 1) then ! attention r may be greater than lo2 or la2 h1 = tab(nx , ny ) h2 = tab(nx +1, ny ) h3 = tab(nx +1, ny +1) h4 = tab(nx , ny +1) hh = h1 * xm * ym + & ! h2 * xp * ym + & ! h3 * xp * yp + & ! h4 * xm * yp ! profile(p + ll / 2 + 1) = hh endif enddo return endsubroutine profile_at_angle endsubroutine plt__acf