Testing posting on Discourse versus Google newsgroup

1

Just comparing posting something from Google newsgroup and Discourse with
something pretty…
snap_59694375-3d42-40d1-ab16-7684fd579430

On a Linux box with fpm, ps2pdf, and evince I ran it with

 fpm run|ps2pdf - >view.pdf;evince view.pdf
see Fortran code here 
xxprogram root3
integer,parameter       :: dp=kind(0.0d0)
real(kind=dp),parameter :: const=1.905e-3_dp, r2r=0.4e0_dp
integer,parameter       :: rows=900, cols=900, r=150
real(kind=dp)           :: angle
complex(kind=dp)        :: z, w, d
integer                 :: i, ix, iy, linebreak
character(len=2)        :: zup
character(len=60)       :: ps(7)=[character(len=60) :: &
&"%!", &
& "/str 256 string def", &
& "/go { image showpage } def", &
& "0 720 translate 72 72 scale", &
& "(I0, 1x, I0, ' 8 [', I0, ' 0 0 ', I0, ' 0 ', I0, ']')", &
& "{currentfile str readhexstring pop}", &
& "go" ]
   write(*,'(A)') (trim(ps(i)),i=1,4)
   write(*,ps(5)) rows, cols, r, -r, r
   write(*,'(A)')(trim(ps(i)),i=6,7)
   linebreak=0
   do ix=-rows/2, rows-rows/2-1
      do iy=-cols/2, cols-cols/2-1
         ! test for cmplx(0.0,0.0) which is undefined
         z=cmplx(ix*const, iy*const, kind=dp)
         do i=1,224
            if ( z == (0.0_dp,0.0_dp) ) exit
            w=z*z
            z=z*w*2+1
            z=z/w/3
            d=z-1
            if(abs(d)<r2r) exit
            d=z+cmplx(0.5_dp, 0.866_dp, kind=dp)
            if(abs(d)<r2r) exit
            d=z+cmplx(0.5_dp, -0.866_dp, kind=dp)
            if(abs(d)<r2r) exit
         enddo
         if(i>224) then
            angle=atan2(0.0_dp,-1.0_dp)
         else
            angle=atan2(d%im, d%re)
         endif
         write(zup,'(Z2.2)') int(abs(angle/0.012272_dp))
         write(*,'(A2)',advance='no') lower(zup)
         linebreak=linebreak+1
         if(linebreak==30)then
            linebreak=0
            write(*,*)
         endif
      enddo
   enddo
contains
elemental pure function lower(str) result (string)
! ident="@(#) M_strings lower(3f) Changes a string to lowercase over specified range"
character(*), intent(in)     :: str
character(len(str))          :: string
integer                      :: i
integer,parameter            :: diff = iachar('A')-iachar('a')
   string = str
   do concurrent (i = 1:len_trim(str))                   ! step thru each letter in the string in specified range
      select case (string(i:i))
      case ('A':'Z')
         string(i:i) = achar(iachar(str(i:i))-diff)   ! change letter to miniscule
      case default
      end select
   enddo
end function lower
end program root3
1 Like
2

Did you see my comments in comp.lang.fortran about this statement? This is just a complicated way to write

angle = 0.0_dp

I suspect that the arguments are switched in the atan2() function. Different languages often have different conventions for this function.

Also this statement

angle=atan2(aimag(d), real(d,kind=dp))

can be written simply as

angle=atan2(d%im, d%re)

There needs to be a test for z=(0.0_dp,0.0_dp) in the loop. Otherwise, the expressions are undefined in fortran and produce NaN on IEEE machines.

Finally, in this version of the code there is

string(i:i) = char(iachar(str(i:i))-diff)

That char() should probably be achar(), right?

3

My port to FORTRAN 77, using SDL77:

sdl77
sdl77814×635 360 KB

FORTRAN 77 
      PROGRAM ROOT3
      EXTERNAL PLOT
      EXTERNAL GCLOSE, GDELAY, GEVENT, GFLUSH, GOPEN
      INTEGER  GKEY

      INTEGER IDELAY, IESC, IEQUIT, IW, IH
      PARAMETER (IDELAY=100, IESC=27, IEQUIT=12, IW=800, IH=600)

      INTEGER IEVENT, ISTAT
      LOGICAL DONE
      DATA DONE /.FALSE./
C
C     OPEN SDL 1.2 WINDOW.
C
      CALL GOPEN(IW, IH, 'FORTRAN' // CHAR(0), ISTAT)
      IF (ISTAT .NE. 0) STOP
C
C     DRAW TO SCREEN LAYER.
C
      CALL GLAYER(0)
      CALL PLOT(IW, IH)
C
C     MAIN LOOP.
C
   10 CONTINUE
   20 CONTINUE
      CALL GEVENT(IEVENT, ISTAT)
      IF (IEVENT .EQ. IEQUIT) DONE = .TRUE.
      IF (ISTAT .EQ. 1) GOTO 20
      IF (GKEY(IESC) .EQ. 1) DONE = .TRUE.
      CALL GFLUSH()
      CALL GDELAY(IDELAY)
      IF (.NOT. DONE) GOTO 10
C
C     CLEAN UP AND QUIT.
C
      CALL GCLOSE()
      END
C     ******************************************************************
      SUBROUTINE PLOT(COLS, ROWS)
      EXTERNAL GCOLOR, GLOCK, GPIXEL, GULOCK

      DOUBLE PRECISION CONST, R2R
      PARAMETER (CONST=1.905D-3, R2R=0.4D0)

      INTEGER COLS, ROWS

      DOUBLE PRECISION ANGLE
      DOUBLE COMPLEX   D, W, Z
      INTEGER          I, IX, IY, J

      CALL GLOCK()

      DO 10 IX = -COLS/2, COLS - COLS/2 - 1
      DO 20 IY = -ROWS/2, ROWS - ROWS/2 - 1
      ANGLE = 0D0
      IF (IX .EQ. 0 .AND. IY .EQ. 0) GOTO 50
      Z = DCMPLX(IX * CONST, IY * CONST)
      DO 30 I = 1, 224
      W = Z**2
      Z = Z * W * 2 + 1
      Z = Z / W / 3
      D = Z - 1
      IF (ABS(D) .LT. R2R) GOTO 40
      D = Z + DCMPLX(0.5D0, 0.866D0)
      IF (ABS(D) .LT. R2R) GOTO 40
      D = Z + DCMPLX(0.5D0, -0.866D0)
      IF (ABS(D) .LT. R2R) GOTO 40
   30 CONTINUE
   40 CONTINUE
      ANGLE = ATAN2(DIMAG(D), DBLE(D))
   50 CONTINUE
      J = INT(ABS(ANGLE / 0.012272D0))
      CALL GCOLOR(J, J, J)
      CALL GPIXEL(IX + COLS/2, IY + ROWS/2)
   20 CONTINUE
   10 CONTINUE

      CALL GULOCK()
      END
3 Likes
4

I intended more to compare posting something graphics-related with example code on the two forums and noted on the comp.lang.fortran I was starting with a verbatim copy of the original post with a few tweeks. The venerable comp.lang.fortran forum has hosted many great discussions and continues to do so. However, for topics such as this, Discourse allows displaying the graphics, highlighted display of the code with proper indenting; and i can correct the code and include those changes. That is appealing and a considerable plus for Discourse when the discussion is not something easily handled from a more than basic terminal-like interface no matter how appealing the simplicity of a comp.lang interface can be (and the Google interface is not as nice as several of the real terminal-based interfaces are, unfortunately. I do grant it is great it is freely available, and feeds are becoming hard to come by, etc.).

Without the original description of the algorithm and having been translated I am sure we can improve it as well, something also done more easily on the now-common modern forum interfaces

The discussion was reminding me of ones begun on the stdlib site about basic support with native Fortran of PPM graphics that has languished. I started a module (M_pixel) along those lines modeled on the VOGLE/M_draw vector graphics library, but needed to add better font support and cleaner line drawing, etc. and really have not proceeded. As noted here, there are some much more capable Fortran-callable interfaces such as SDL77, so I am pondering that as well; as I see pixmap examples in Fortran that always have to resort to PostScript or external apps and think there is probably a place for a basic pixel-based Fortran library, but not sure what form it should take, etc. But for the purposes here; I can update the code in Discourse. Not is not the case in a newsgroup; where you have to repost. Not at a particularly good machine right now, but I will include the changes soon. Thanks! Like the original post, not sure there is any use for it except as an example of the beauty sometimes hidden in simple mathematics, but I like it and will also probably keep a copy, so it is nice to clean it up.

The arguments in C and Fortran for atan2 are both (y,x), although I would prefer (x,y) myself, and have caught myself assuming more than once! Definitely achar().
The atan2 is probably trying to return PI, so that needs changed; order and sign issues there I think. You were right to be suspicious. Thinking about it.
Especially since original C had an unused PI value.

Think I got them @RonShepard . Please take a look. Maybe too obtuse on the cmplx(0.,0.)?

PS: Given in Fortran that

expression            | value
atan2(1.0_dp,0.0_dp)  |  1.5707963267948966
atan2(0.0_dp,1.0_dp)  |  0.0000000000000000
atan2(-1.0_dp,0.0_dp) | -1.5707963267948966
atan2(0.0_dp,-1.0_dp) |  3.1415926535897931
cmplx(0.0_dp,0.0_dp)  |  NAN

those changes make sense to me.

FOOTNOTE

With the suggested changes, the output of the original C program posted at comp.lang.fortran and the program are identical; they were not without the changes. So at least the bugs should now be consistent :slight_smile:

5

I think the test should be inside the do i=1,224 loop. In the code I posted at c.l.f, the test I used was

if ( z == (0.0_dp,0.0_dp) ) exit

The potential problem is three statements later, z=z/w/3.

Regarding atan2() in various languages, in excel the arguments are switched, so I remain constantly confused. In my own fortran codes, I often use keyword arguments for that function, atan2(x=xval,y=yval), so I know I always get it right that way. I presume that whoever first translated the fortran (or C) code might have gotten those arguments switched.

7

Was the convention of y first for atan2(y,x) chosen because if x>0 then atan2(y,x) = atan(y/x) ?

I have seen mathematics textbooks that say the principal argument of the complex number x+iy is atan(y/x), which is bad if x=0 and false if x<0 with the usual definition of “principal”. When teaching complex analysis I tried telling the students that branch cuts are like the International Date Line: the branch points are fixed by the function one is dealing with but the path the cut follows between them depends on what problem is being dealt with. For the Date Line the function is local time, the branch points are the poles, and the problem is politics. Unfortunately that didn’t help the many students who hadn’t understood the Date Line when they heard about it at school and hadn’t crossed it :frowning:

1 Like