%pl issue

The program below uses %pl to plot the Bessel function Y1(x) in the range 0 < x < 20

Obviously the first point in the x data arrays cannot be zero since Y1(0) is negative infinity.

We can also use the %pl option y_min to clip the large negative values (so effectively the line being plotted starts somewhere outside the %pl frame, and should appear where it crosses the lower horizontal edge of the frame.

Certain combinations of the starting x value (and hence starting y value) combined with the value specified for y_min run into problems. This can be seen in the second plot produced by this program. Here is a screenshot.

https://www.dropbox.com/scl/fi/fodkuemt1ykexr10v2vl2/Screenshot-2026-01-09-105813.jpg?rlkey=3u7gogu8lbem1xw52ejkbu9re&st=4ld4pjgp&dl=0

program test
use clrwin
implicit none
real*8 x_pl1(200), y1_pl1(200)
real*8 x_pl2(200), y1_pl2(200)
real*8 x_pl3(200), y1_pl3(200)
  integer i, iw
    x_pl1(1) = 0.0001d0   !Y0 of zero is - infinity so start a small distance away from zero
    do i = 2, 200
      x_pl1(i) = x_pl1(i-1) + 0.1d0
    end do
    y1_pl1=bessel_y1(x_pl1)
    print*, 'First point in first graph at : ', x_pl1(1), y1_pl1(1)
    call winop@('%pl[native,n_graphs=1,x_array,width=3]')
    call winop@('%pl[frame,etched,gridlines]')
    call winop@('%pl[link=lines, colour=blue]')
    call winop@('%pl[y_max=1,y_min=-1,x_max=20]')   ! Clip y-axis at -1
    iw = winio@('%pl&',500,400,200,x_pl1,y1_pl1)
    
    x_pl2(1) = 0.00001d0 ! Start closer to zero
    do i = 2, 200
      x_pl2(i) = x_pl2(i-1) + 0.1d0
    end do
    y1_pl2=bessel_y1(x_pl2)
    print*, 'First point in second graph at : ', x_pl2(1), y1_pl2(1)
    call winop@('%pl[native,n_graphs=1,x_array,width=3]')
    call winop@('%pl[frame,etched,gridlines]')
    call winop@('%pl[link=lines,colour=blue]')
    call winop@('%pl[y_max=1,y_min=-1,x_max=20]') ! Clip y-axis at -1
    iw = winio@('%pl&',500,400,200,x_pl2,y1_pl2)

    x_pl3(1) = 0.00001d0
    do i = 2, 200
      x_pl3(i) = x_pl3(i-1) + 0.1d0
    end do
    y1_pl3=bessel_y1(x_pl3)
    print*, 'First point in third graph at : ', x_pl3(1), y1_pl3(1)
    call winop@('%pl[native,n_graphs=1,x_array,width=3]')
    call winop@('%pl[frame,etched,gridlines]')
    call winop@('%pl[link=lines,colour=blue]')
    call winop@('%pl[y_max=1,y_min=-2,x_max=20]') ! Clip y-axis at -2
    iw = winio@('%pl',500,400,200,x_pl3,y1_pl3)
end program test

PS If i plot -Y1, i.e. the mirror image (reflected in the X axis) everthing is ok.

Thank you for this feedback. I will take a look at it.

Had a quick look and setting these values the same, removes the issue with the second plot:

x_pl1(1) = 0.0001d0 x_pl2(1) = 0.0001d0 x_pl3(1) = 0.0001d0

First point in first graph at : 1.000000000000E-04 -6366.19803646 First point in second graph at : 1.000000000000E-04 -6366.19803646 First point in third graph at : 1.000000000000E-04 -6366.19803646

Does this look correct?

Lester

Lester, the point Kenneth was making was 'It shouldn't have done it like this.'

I think we all have had moments like this, where what we did seems to violate some unknown limitation, so we programmed around it. Perfectly fine, yet the issue still remains.

Like my post a few days ago about property sheets locking up and crashing. Shouldn't happen, yet it does. So we do something different.

The folks at Silverfrost do a great job with this product, and knowing that there is something 'odd' is sometimes enough to get an even more robust product. If I dug back in my memory, I could give you at LEAST a half dozen workarounds in both MS Visual C++ and in Silverfrost. Only Silverfrost responded positively to my bug reports. And I learned a bunch of new things as a result of many of these encounters.

I hope that this is just one of those 'oh, goodness, look at this!' moments where things get corrected and more robust.

I wholeheartedly agree with Blll’s sentiments “Silverfrost do a great job with this product”.

I hope new readers of the forum don’t misinterpret my many posts on issues with %pl i.e. getting the impression that it is a “fragile” element of Clearwin. I have been pushing %pl far beyond a simple x-y plotting procedure. There are many advanced features available to the Fortran programmer using Clearwin and %pl. For example, after reading the Wikipedia article on domain colouring on Saturday morning, by early Saturday evening I had a working domain_colour_plot. This experimental code below demonstrates how powerful %pl really is!

module domain_plot_mod
use clrwin
implicit none
private
public domain_plot

real*8, allocatable:: xc(:), yc(:), amp(:,:), phase(:,:)
real*8, parameter :: pi = 4.d0*atan(1.d0)
real*8, parameter :: twopi = 2.d0*pi
real*8, parameter :: halfpi = 0.5d0*pi
real*8, parameter :: rec2pi = 1.d0/(2.d0*pi)
real*8 :: xmin, xmax, ymin, ymax, sat, zabsmax

contains

  subroutine domain_plot(x,y,z,s)
  !------------------------------------------------------------------
  ! domain_plot plots a complex function f(z) over a 2D domain,
  ! mapping magnitude to brightness and phase to hue, with automatic
  ! magnitude and phase contour highlighting. Contour widths decrease
  ! in proportion to the local rate of change.
  !------------------------------------------------------------------
  real*8, intent(in) :: x(:)
  real*8, intent(in) :: y(:)
  complex*16, intent(in) :: z(:,:)
  real*8 :: s
  integer :: i, j, nx, ny, iw
    nx = size(x) ; ny = size(y)
    if (size(z,dim=1) .ne. nx .or. size(z,dim=2) .ne. ny) &
      STOP 'Non-conformant arrays in call to DOMAIN_PLOT'
    allocate (xc(nx),yc(nx),amp(nx,ny),phase(nx,ny))
    xmin = minval(x) ; xmax = maxval(x)
    ymin = minval(y) ; ymax = maxval(y)
    xc = x ; yc = y ; sat = s
    sat = max(0.d0,(min(sat,1.d0)))
    zabsmax = 0.d0
    do i = 1, nx
      do j = 1, ny
        amp(i,j)   = abs(z(i,j))
        zabsmax = max(zabsmax,amp(i,j))
        if (amp(i,j) > 0.d0) then
          phase(i,j) = atan2(aimag(z(i,j)), real(z(i,j)))
        else
          phase(i,j) = 0.d0
        end if
      end do
    end do
    if (zabsmax > 0.d0) then
      iw = winio@('%fn[Arial]%ts%bf%bg&',1.2d0,rgb@(220,220,220))
      call winop@('%pl[native, n_graphs=1,x_array,frame,link=user,external_ticks]')
      call winop@('%pl[axes_pen=2,frame_pen=2,xaxis="Re",yaxis="Im"]')
      call winop_flt@('%pl[x_min]',xmin)
      call winop_flt@('%pl[x_max]',xmax)
      call winop_flt@('%pl[y_min]',ymin)
      call winop_flt@('%pl[y_max]',ymax) 
      iw = winio@('%^pl',800,600,2,[0.d0,1.d0],[0.d0,1.d0],pl_cb)
    else
      print*, 'Domain_plot: Input is level plane.'
    end if
    deallocate(xc,yc,amp,phase)
  end subroutine domain_plot

  integer function pl_cb()
  integer :: ixmin, ixmax, iymin, iymax, i, j, k, icol
  real*8 :: rxmin, rxmax, rymin, rymax, x, y, r, h
  real*8 :: red, green, blue, t, f, d, line
  real*8 :: target_phase, fphase
  real*8 :: v, scale_factor
  real*8, parameter :: log2 = log(2.d0)
  real*8, parameter :: inv_log2 = 1.d0/log2
  real*8, parameter :: width = 0.06d0       ! a fraction of a log-magnitude decade
  real*8, parameter :: strength = 0.35d0
  real*8, parameter :: ref_mag = 1000.d0
  real*8, parameter :: v_min = 0.25d0       ! minimum brightness to prevent black pixels
    scale_factor = ref_mag/zabsmax
    if (clearwin_string@('CALLBACK_REASON') .eq. 'PLOT_ADJUST') then
      k = GET_PLOT_POINT@(xmin,ymin,rxmin,rymin)
      k = GET_PLOT_POINT@(xmax,ymax,rxmax,rymax)
      ixmin = nint(rxmin) ; ixmax = nint(rxmax)
      iymin = nint(rymin) ; iymax = nint(rymax)
      do i = ixmin, ixmax, 1
        do j = iymin, iymax, -1
          k = get_plot_data@(i,j,x,y)
          ! Interpolate amplitude and phase
          call z_rh_bilinear(xc,yc,amp,phase,x,y,r,h)
          h = (h + pi) * rec2pi
          h = h - floor(h)
          r = r*scale_factor
          ! Base brightness from amplitude
          v = 1.d0 - 1.d0/(1.d0 + log(1.d0 + r))
          v = max(v, v_min)
          ! Magnitude contours (white rings)
          if (r > 0.d0) then
            t = log(r)*inv_log2
            f = t - floor(t)
            d = min(f, 1.d0 - f)
            line = exp(-(d/width)**2)
            ! Blend towards white
            v = v + (1.d0 - v)*line*strength
          end if
          ! Phase contours (white lines)
          do k = 0, 3
            target_phase = halfpi*dble(k)
            fphase = abs(h*twopi - target_phase)
            if (fphase > pi) fphase = twopi - fphase
            line = exp(-(fphase/width)**2)
            ! Blend towards white
            v = v + (1.d0 - v)*line*strength
          end do
          ! Final clamp of brightness
          v = min(1.d0, max(v_min, v))
          ! Convert HSV to RGB
          call hsv_to_rgb(h,sat,v,red,green,blue)
          icol = rgb@(int(255.d0*red), int(255.d0*green), int(255.d0*blue))
          ! Draw pixels using draw_filled_rectangle@
          call draw_filled_rectangle@(i,j,i,j,icol)
        end do
      end do
    end if
    pl_cb = 2
  end function pl_cb

  subroutine z_rh_bilinear(x, y, amp, phase, xt, yt, r, h)
  real*8, intent(in)  :: x(:), y(:)
  real*8, intent(in)  :: amp(:,:), phase(:,:)
  real*8, intent(in)  :: xt, yt
  real*8, intent(out) :: r, h
  integer :: i, j
  real*8 :: x1, x2, y1, y2, tx, ty
  real*8 :: a11, a21, a12, a22
  real*8 :: c11, c21, c12, c22
  real*8 :: s11, s21, s12, s22
  real*8 :: c, s, rrr
    ! Locate cell
    i = bin_search(x, xt) ; j = bin_search(y, yt)
    x1 = x(i)   ;  x2 = x(i+1)
    y1 = y(j)   ;  y2 = y(j+1)
    tx = (xt - x1) / (x2 - x1)
    ty = (yt - y1) / (y2 - y1)
    ! Amplitude at corners
    a11 = amp(i,  j  )
    a21 = amp(i+1,j  )
    a12 = amp(i,  j+1)
    a22 = amp(i+1,j+1)
    ! Phase as unit complex numbers
    c11 = cos(phase(i,  j  )) ; s11 = sin(phase(i,  j  ))
    c21 = cos(phase(i+1,j  )) ; s21 = sin(phase(i+1,j  ))
    c12 = cos(phase(i,  j+1)) ; s12 = sin(phase(i,  j+1))
    c22 = cos(phase(i+1,j+1)) ; s22 = sin(phase(i+1,j+1))
    ! Bilinear interpolation of amplitude
    r = (1.d0-tx)  *(1.d0-ty)*a11 + &
         tx        *(1.d0-ty)*a21 + &
          (1.d0-tx)*ty       *a12 + &
          tx       *ty       *a22
    ! Bilinear interpolation on unit circle
    c = (1.d0-tx)  *(1.d0-ty)*c11 + &
         tx        *(1.d0-ty)*c21 + &
         (1.d0-tx) *ty       *c12 + &
         tx        *ty       *c22
    s = (1.d0-tx)  *(1.d0-ty)*s11 + &
         tx        *(1.d0-ty)*s21 + &
         (1.d0-tx) *ty        *s12 + &
         tx        *ty        *s22
    ! Normalise and recover phase
    rrr = sqrt(c*c + s*s)
    if (rrr > 0.d0) then
      h = atan2(s/rrr, c/rrr)
    else
      h = 0.d0
    end if
  end subroutine z_rh_bilinear

  integer function bin_search(vec, val)
  real*8, intent(in) :: vec(:), val
  integer :: lo, hi, mid
    lo = 1 ; hi = size(vec)-1
    do while (lo .le. hi)
      mid = (lo + hi) / 2
      if (val .lt. vec(mid)) then
        hi = mid - 1
      else if (val .gt. vec(mid+1)) then
        lo = mid + 1
      else
        bin_search = mid
        return
      end if
    end do
    ! If outside range, clamp to last valid index
    bin_search = max(1, min(size(vec)-1, lo))
  end function bin_search

  subroutine hsv_to_rgb(h, s, v, r, g, b)
  real*8, intent(in) :: h, s, v
  real*8, intent(out) :: r, g, b
  real*8 :: i, f, p, q, t
    i = floor(6.d0*h)
    f = 6.d0*h - i
    p = v*(1.d0 - s)
    q = v*(1.d0 - f*s)
    t = v*(1.d0 - (1.d0-f)*s)
    select case (mod(int(i),6))
      case (0); r=v; g=t; b=p
      case (1); r=q; g=v; b=p
      case (2); r=p; g=v; b=t
      case (3); r=p; g=q; b=v
      case (4); r=t; g=p; b=v
      case (5); r=v; g=p; b=q
    end select
  end subroutine hsv_to_rgb

end module domain_plot_mod

program test_domain_plot
use domain_plot_mod, only : domain_plot
implicit none
integer, parameter :: dp = kind(1.d0)
integer, parameter :: n = 200
integer :: i, j
real(dp) :: re_z(n), im_z(n)
complex(dp) :: f(n,n), z_temp, cj
  ! Generate first test function
  cj = cmplx(0.d0,1.d0,kind=dp)
  do i = 1, n
    re_z(i) = -2.0d0 + 5.0d0 * dble(i-1) / dble(n-1)
    do j = 1, n
      im_z(j)  = -2.0d0 + 5.0d0 * dble(j-1) / dble(n-1)
      z_temp = cmplx(re_z(i),im_z(j),kind=dp)
      f(i,j) = ((z_temp*z_temp - 1.d0)*(z_temp - 2.d0 - cj)**2 )/(z_temp*z_temp+2.d0+2.d0*cj)
    end do
  end do
  ! Plot function with different saturation levels applied
  call domain_plot(re_z,im_z,f,0.d0)
  call domain_plot(re_z,im_z,f,0.5d0)
  call domain_plot(re_z,im_z,f,1.0d0)
  ! Generate second test functon
  do i = 1, n
    re_z(i) = -3.d0 + 6.d0*dble(i-1)/dble(n-1)
    do j = 1, n
      im_z(j) = -3.d0 + 6.d0*dble(j-1)/dble(n-1)
      z_temp = cmplx(re_z(i),im_z(j),kind=dp)
      f(i,j) = (z_temp-1.d0)/((z_temp**2)+z_temp+1.d0)
    end do
  end do
  ! Plot the second test function with different saturation levels applied
  call domain_plot(re_z,im_z,f,0.d0)
  call domain_plot(re_z,im_z,f,0.25d0)
  call domain_plot(re_z,im_z,f,0.80d0)
end program test_domain_plot

This is a truly impressive piece of programming given the time line for development.

The runtime needed to create the images is naturally quite long because of the pixel-by-pixel rendering. If this is an issue of concern then one could consider trying a) a ClearWin+ "user_surface" and/or b) some kind of parallel processing.

Thanks Paul, Now much improved and in the sample code sub-forum.
I better think about doing some real work next week!