I have always been taught that blue is the color of lies because you
have fewer cones sensitive to blue so that you don't get as clear a
picture in blue as in red or green. So here is a Fortran program that
permits testing. It spits out 4 images, in red, green, blue, and gray
in *.PPM format, which is about the simplest image format possible.
Unfortunately, the only software I have that can open *.PPM files
is gimp, which is free. Gimp can also export to other, more common
formats.

program lies
use ISO_FORTRAN_ENV
implicit none
integer iunit
integer i, j, k
character(9):: filename(4)=[character(9)::'red.PPM','green.PPM', &
'blue.PPM','gray.PPM']
integer:: color(3,4)=reshape([1,0,0,0,1,0,0,0,1,1,1,1],[3,4])
do i = 1, 4
open(newunit=iunit,file=filename(i),status='replace', &
access='stream')
write(iunit) 'P6 3840 2160 255 '
write(iunit) ((int((j-1)/15*color(:,i),INT8),j=1,3840),k=1,2160)
close(iunit)
end do
end program lies

In which color is it easiest to identify the changes in intensity? In
which is it the most difficult, i.e., the color of lies?

Not really all that Fortran-related, but I was surprised how brief
a program can be that creates these UHD images.

On Wednesday, August 31, 2022 at 12:07:49 AM UTC-7, James Van Buskirk wrote:
> I have always been taught that blue is the color of lies because you
> have fewer cones sensitive to blue so that you don't get as clear a
> picture in blue as in red or green.

It isn't Fortran, but I wrote this one some years ago. It generates PostScript
directly, easy to print.

For each pixel in the output image, it computes a complex cube root.

I have thought about doing the whole calculation in PS, but didn't do that yet.

#include <stdio.h>
#include <math.h>

#define MULT(x,y,a,b,t) ((t)=(x)*(a)-(y)*(b),(y)=(x)*(b)+(y)*(a),(x)=(t))
#define DIV(x,y,a,b,t) (MULT(x,y,a,-b,t),x /= SQ(a,b),y /= SQ(a,b))
#define SQ(x,y) ((x)*(x)+(y)*(y))
#define pi 3.141592654
#define R 150

#define CONST 1.905e-3
#define ROWS 900
#define COLS 900
#define R2 .16

main() {
double x,y,zr,zi,wr,wi,dr,di,angle,t;
int i,j,gray,ix,iy;
static char dots[]={0x80,0x40,0x20,0x10,0x08,0x04,0x02,0x01};
static char *ps[]={"%%!\n",
"/str 256 string def\n",
"/go { image showpage } def\n",
"0 720 translate 72 72 scale\n",
"%d %d 8 [%d 0 0 -%d 0 %d]\n",
"{currentfile str readhexstring pop}\n",
"go\n"};

for(i=0;i<sizeof(ps)/sizeof(char*);i++) printf(ps[i],ROWS,COLS,R,R,R);

j=0;
for(ix=-ROWS/2;ix<ROWS-ROWS/2;ix++)
for(iy=-COLS/2;iy<COLS-COLS/2;iy++) {
zr=ix*CONST;
zi=iy*CONST;
gray=65535;
for(i=1;i<=224;i++) {
wr=zr,wi=zi;
MULT(wr,wi,zr,zi,t);
MULT(zr,zi,wr,wi,t);
zr *=2, zi *= 2;
zr += 1;
DIV(zr,zi,wr,wi,t);
zr /= 3; zi /= 3;
dr=zr-1,di=zi;
if(SQ(dr,di)<R2) break;
dr=zr+0.5,di=zi+0.866;
if(SQ(dr,di)<R2) break;
dr=zr+0.5,di=zi-0.866;
if(SQ(dr,di)<R2) break;
}
if(i>224) angle=atan2(0.,-1.);
else angle=atan2(di,dr);
#if 0
fprintf(stderr,"%d %d %f %d\n",ix,iy,angle,i);
#endif
printf("%02x",(int)(fabs(angle)/.012272));
if(++j%30==0) putchar('\n');
}
}

Well, this isn't Fortran either, but the Tcl program below displays part of the PPM files (my screen is too small to show it all):

# viewppm.tcl --
# View a PPM file
# pack [canvas .c -bg white]

set image [image create photo -file green.PPM]

..c create image 0 0 -image $image -anchor c

I use an anchor "central" because at first I only saw the top-left, which looks like a uniform black rectangle.

Regards,

Arjen

On Wednesday, August 31, 2022 at 12:07:49 AM UTC-7, James Van Buskirk wrote:
> I have always been taught that blue is the color of lies because you
> have fewer cones sensitive to blue so that you don't get as clear a
> picture in blue as in red or green.

OK, here it is in Fortran. The original was before C had complex,
and it looks much nicer with complex variables.
(Note that it is fixed form, and many news readers lose leading blanks.)

PROGRAM ROOT3
REAL*8 PI/3.141592653589793/, CONST/1.905E-3/, R2R/.4/
INTEGER ROWS/900/, COLS/900/, R/150/
REAL*8 ANGLE
COMPLEX*16 Z, W, D
INTEGER I, J, GRAY, IX, IY
CHARACTER*50 PS(7)
DATA PS/"%!",
* "/str 256 string def",
* "/go { image showpage } def",
* "0 720 translate 72 72 scale",
* "(I6, I6, ' 8 [', I6, ' 0 0 ', I6, ' 0 ', I6, ']')",
* "{currentfile str readhexstring pop}",
* "go"/
DO I=1,4
WRITE(*,'(A)') PS(I)
END DO
WRITE(*,PS(5)) ROWS, COLS, R, -R, R
DO I=6,7
WRITE(*,'(A)') PS(I)
END DO

DO IX=-ROWS/2, ROWS-ROWS/2-1
DO IY=-COLS/2, COLS-COLS/2-1
Z=CMPLX(IX*CONST, IY*CONST, KIND(1.D0))
DO I=1,224
W=Z*Z
Z=Z*W*2+1
Z=Z/W/3
D=Z-1
IF(ABS(D)<R2R) EXIT
D=Z+CMPLX(0.5D0, 0.866D0, KIND(1.D0))
IF(ABS(D)<R2R) EXIT
D=Z+CMPLX(0.5D0, -0.866D0, KIND(1.D0))
IF(ABS(D)<R2R) EXIT
END DO
IF(I>224) THEN
ANGLE=ATAN2(0.,1.)
ELSE
ANGLE=ATAN2(AIMAG(D), REAL(D))
ENDIF
WRITE(*,'(Z2)') INT(ABS(ANGLE/0.012272))
END DO
END DO
END

The beauty that can emerge from math can be startling. Nice image.
It was fun playing with the formula and seeing the different results.
Just a note that to work in my programming environment I had to change
the Z2 format to Z2.2.

On Wednesday, August 31, 2022 at 7:59:12 PM UTC-7, John wrote:
> The beauty that can emerge from math can be startling. Nice image.
> It was fun playing with the formula and seeing the different results.
> Just a note that to work in my programming environment I had to change
> the Z2 format to Z2.2.

The C version writes 60 hex digits per line, but it isn't so easy to do in Fortran.
Could be done with non-advancing I/O, but I did try that. PostScript doesn't care.

And yes, it does need Z2.2. I was testing different things before posting,
and didn't get that one. I am used to %02x, but not so used to that one.

On Thursday, September 1, 2022 at 8:12:22 AM UTC+10, gah4 wrote:
> On Wednesday, August 31, 2022 at 12:07:49 AM UTC-7, James Van Buskirk wrote:
> > I have always been taught that blue is the color of lies because you
> > have fewer cones sensitive to blue so that you don't get as clear a
> > picture in blue as in red or green.
> OK, here it is in Fortran. The original was before C had complex,
> and it looks much nicer with complex variables.
> (Note that it is fixed form, and many news readers lose leading blanks.)
>
> PROGRAM ROOT3
> REAL*8 PI/3.141592653589793/, CONST/1.905E-3/, R2R/.4/

Why is it that experienced programmers still get this wrong
after 50 years? All those constants need D0, D-3, D0 respectively.
You should use a decent language -- PL/I.

What's more, REAL*8 is non-standard Fortran.

> INTEGER ROWS/900/, COLS/900/, R/150/
> REAL*8 ANGLE
> COMPLEX*16 Z, W, D
> INTEGER I, J, GRAY, IX, IY
> CHARACTER*50 PS(7)
> DATA PS/"%!",
> * "/str 256 string def",
> * "/go { image showpage } def",
> * "0 720 translate 72 72 scale",
> * "(I6, I6, ' 8 [', I6, ' 0 0 ', I6, ' 0 ', I6, ']')",
> * "{currentfile str readhexstring pop}",
> * "go"/
> DO I=1,4
> WRITE(*,'(A)') PS(I)
> END DO
> WRITE(*,PS(5)) ROWS, COLS, R, -R, R
> DO I=6,7
> WRITE(*,'(A)') PS(I)
> END DO
>
> DO IX=-ROWS/2, ROWS-ROWS/2-1
> DO IY=-COLS/2, COLS-COLS/2-1
> Z=CMPLX(IX*CONST, IY*CONST, KIND(1.D0))
> DO I=1,224
> W=Z*Z
> Z=Z*W*2+1
> Z=Z/W/3
> D=Z-1
> IF(ABS(D)<R2R) EXIT
> D=Z+CMPLX(0.5D0, 0.866D0, KIND(1.D0))
> IF(ABS(D)<R2R) EXIT
> D=Z+CMPLX(0.5D0, -0.866D0, KIND(1.D0))
> IF(ABS(D)<R2R) EXIT
> END DO
> IF(I>224) THEN
> ANGLE=ATAN2(0.,1.)
> ELSE
> ANGLE=ATAN2(AIMAG(D), REAL(D))
> ENDIF
> WRITE(*,'(Z2)') INT(ABS(ANGLE/0.012272))
> END DO
> END DO
> END

On Thursday, September 1, 2022 at 12:12:22 AM UTC+2, gah4 wrote:
> On Wednesday, August 31, 2022 at 12:07:49 AM UTC-7, James Van Buskirk wrote:
> > I have always been taught that blue is the color of lies because you
> > have fewer cones sensitive to blue so that you don't get as clear a
> > picture in blue as in red or green.
> OK, here it is in Fortran. The original was before C had complex,
> and it looks much nicer with complex variables.
> (Note that it is fixed form, and many news readers lose leading blanks.)
>

I tried to view the result of this Fortran program with ghostscript and got rubbish. I also tried the result of the C program and that showed the picture I expected. So there is something different, possibly caused by the formatting. But I am no expert in PostScript.

Regards,

Arjen

On Wednesday, August 31, 2022 at 11:49:35 PM UTC-7, Robin Vowels wrote:

(snip)

> > REAL*8 PI/3.141592653589793/, CONST/1.905E-3/, R2R/.4/
> Why is it that experienced programmers still get this wrong
> after 50 years? All those constants need D0, D-3, D0 respectively.

They really don't, and even more probably don't need to be REAL*8.

Especially as PI isn't even used in the program!

In the end, the program converts an angle to a gray level from 0 to 255,
which the printer likely converts to a coarser gray level.

CONST relates to the scale of the image, which also doesn't need
to be very exact.

On Thursday, September 1, 2022 at 1:39:32 AM UTC-7, arjen.m...@gmail.com wrote:

(snip)

> I tried to view the result of this Fortran program with ghostscript and got rubbish.
> I also tried the result of the C program and that showed the picture I expected.
> So there is something different, possibly caused by the formatting.
> But I am no expert in PostScript.

Yes, as John noted the format should be Z2.2.

It is not so obvious, but PS reads the hex characters ignoring any other characters.

With Z2 it blank pads, so only one digit. Z2.2 zero pads, so it works.
The C program uses %02x, which also zero pads.

On Thursday, September 1, 2022 at 6:45:50 PM UTC+10, gah4 wrote:
> On Wednesday, August 31, 2022 at 11:49:35 PM UTC-7, Robin Vowels wrote:
>
> (snip)
> > > REAL*8 PI/3.141592653589793/, CONST/1.905E-3/, R2R/.4/
> > Why is it that experienced programmers still get this wrong
> > after 50 years? All those constants need D0, D-3, D0 respectively.
> They really don't, and even more probably don't need to be REAL*8.
..
It's sloppy programming, and a bad example for others.
..
> Especially as PI isn't even used in the program!
>
> In the end, the program converts an angle to a gray level from 0 to 255,
> which the printer likely converts to a coarser gray level.
>
> CONST relates to the scale of the image, which also doesn't need
> to be very exact.

On 9/1/2022 3:57 AM, gah4 wrote:
> On Thursday, September 1, 2022 at 1:39:32 AM UTC-7, arjen.m...@gmail.com wrote:
>
> (snip)
>
>> I tried to view the result of this Fortran program with ghostscript and got rubbish.
>> I also tried the result of the C program and that showed the picture I expected.
>> So there is something different, possibly caused by the formatting.
>> But I am no expert in PostScript.
>
> Yes, as John noted the format should be Z2.2.
>
> It is not so obvious, but PS reads the hex characters ignoring any other characters.
>
> With Z2 it blank pads, so only one digit. Z2.2 zero pads, so it works.
> The C program uses %02x, which also zero pads.
>
>

Another tiny correction to the program:

When ix = iy = 0, the statement

Z=Z/W/3

causes division of complex zero by complex zero to occur. With some
Fortran compilers, this might result in '**' to be printed out
instead of two hex digits; in that case, a Postscript interpreter
may abort the display of the image file.

-- CyrMag

On Thursday, September 1, 2022 at 4:31:55 AM UTC-7, CyrMag wrote:
> On 9/1/2022 3:57 AM, gah4 wrote:

(snip)

> > Yes, as John noted the format should be Z2.2.
> > It is not so obvious, but PS reads the hex characters ignoring any other characters.
> > With Z2 it blank pads, so only one digit. Z2.2 zero pads, so it works.
> > The C program uses %02x, which also zero pads.
> Another tiny correction to the program:
> When ix = iy = 0, the statement
> Z=Z/W/3
> causes division of complex zero by complex zero to occur. With some
> Fortran compilers, this might result in '**' to be printed out
> instead of two hex digits; in that case, a Postscript interpreter
> may abort the display of the image file.
I always wondered about the 224 loop, and the test at the end.

It seems that the one case where it happens is the case you mention.
So, angle = atan(1.,0.) = 0, and it works fine.

As well as I know, readhexstring ignores any character that isn't a hex digit.
There would, then, in the case of ** be a one pixel offset for the rest of the
image, which might be noticeable.

Not counting that case, the highest I seems to be 51.

I am not sure now where this is from. I believe I wrote the C program
from a description of the algorithm. I didn't select the 224.

Postscript interpreters I know, if image gives up early, will execute
the showpage and print the image so far. They might not all do that.

Since it came up, hopefully this demonstrates posting in a newsgroup
with indenting ...

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

fpm run|ps2pdf - >view.pdf;evince view.pdf

If indented four or more characters that line should show as-is, but a ">" is
more reliable for verbatim text.

This is the same program essentially verbatim (hopefully) with the
advance='no' added and the Z values lowercased as I was doing a diff on
the output with the C version. So it was tweeked to write thirty values
per line in lowercase to make that easier, and uses the Z2.2 format,
and shows that (hopefully) newgroup viewers leave lines starting with ">"
as-is, which should retain indenting. Everything else is just cosmetic to
satisfy my personal preferences. It does not address the Z/W/3 issue. Like
the original post (with the Z2.2 format) it worked with gfortran, ifort,
and nvfortran.

> program 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
> z=cmplx(ix*const, iy*const, kind=dp)
> do i=1,224
> 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(aimag(d), real(d,kind=dp))
> 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
> character(*), intent(in) :: str
> character(len(str)) :: string
> integer :: i
> integer,parameter :: diff = iachar('A')-iachar('a')
> string = str
> ! step thru each letter in the string in specified range
> do concurrent (i = 1:len_trim(str))
> select case (string(i:i))
> case ('A':'Z')
> string(i:i) = char(iachar(str(i:i))-diff) ! change letter to miniscule
> case default
> end select
> enddo
> end function lower
> end program root3

Well, that worked for a long time, but google groups at least still shows no indenting. Oops. That thing goes out of it's way to mess up any kind of formatting, etc.

On Thursday, September 1, 2022 at 5:05:14 PM UTC-4, John wrote:
> Since it came up, hopefully this demonstrates posting in a newsgroup
> with indenting ...
>
> On a Linux box with fpm, ps2pdf, and evince I ran it with
>
> fpm run|ps2pdf - >view.pdf;evince view.pdf
>
> If indented four or more characters that line should show as-is, but a ">" is
> more reliable for verbatim text.
>
> This is the same program essentially verbatim (hopefully) with the
> advance='no' added and the Z values lowercased as I was doing a diff on
> the output with the C version. So it was tweeked to write thirty values
> per line in lowercase to make that easier, and uses the Z2.2 format,
> and shows that (hopefully) newgroup viewers leave lines starting with ">"
> as-is, which should retain indenting. Everything else is just cosmetic to
> satisfy my personal preferences. It does not address the Z/W/3 issue. Like
> the original post (with the Z2.2 format) it worked with gfortran, ifort,
> and nvfortran.
>
> > program 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
> > z=cmplx(ix*const, iy*const, kind=dp)
> > do i=1,224
> > 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(aimag(d), real(d,kind=dp))
> > 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
> > character(*), intent(in) :: str
> > character(len(str)) :: string
> > integer :: i
> > integer,parameter :: diff = iachar('A')-iachar('a')
> > string = str
> > ! step thru each letter in the string in specified range
> > do concurrent (i = 1:len_trim(str))
> > select case (string(i:i))
> > case ('A':'Z')
> > string(i:i) = char(iachar(str(i:i))-diff) ! change letter to miniscule
> > case default
> > end select
> > enddo
> > end function lower
> > end program root3
So doing anything with graphics I think I am going to go with Discourse

https://fortran-lang.discourse.group/t/testing-posting-on-discourse-versus-google-newsgroup/4256

On Thursday, September 1, 2022 at 10:39:32 AM UTC+2, Arjen Markus wrote:

> >
> I tried to view the result of this Fortran program with ghostscript and got rubbish. I also tried the result of the C program and that showed the picture I expected. So there is something different, possibly caused by the formatting. But I am no expert in PostScript.
>

Since I get a completely contorted picture with the Fortran program I am trying to figure out what is going on. The results are identical between Intel Fortran oneAPI and gfortran, but differ greatly when using the C version.

I have tried analysing the differences between the output from the Fortran program and the C program. So far I can find no cause, but:
- The output file from the C program is roughly 54000 bytes larger
- The number of lines is the same (I used a similar modification as John, but by filling a string of 60 characters)
- The lengths of the lines are equal
- The differences occur in the hexadecimal output but curiously enough not in the actual data (I write out the variables)

Very odd. I must be looking at it in a wrong way.

Regards,

Arjen

On 9/1/2022 11:24 PM, Arjen Markus wrote:
> On Thursday, September 1, 2022 at 10:39:32 AM UTC+2, Arjen Markus wrote:
>
>> I tried to view the result of this Fortran program with ghostscript and got rubbish. I also tried the result of the C program and that showed the picture I expected. So there is something different, possibly caused by the formatting. But I am no expert in PostScript.
>>
> Since I get a completely contorted picture with the Fortran program I am trying to figure out what is going on. The results are identical between Intel Fortran oneAPI and gfortran, but differ greatly when using the C version.
>
> I have tried analysing the differences between the output from the Fortran program and the C program. So far I can find no cause, but:
> - The output file from the C program is roughly 54000 bytes larger
> - The number of lines is the same (I used a similar modification as John, but by filling a string of 60 characters)
> - The lengths of the lines are equal
> - The differences occur in the hexadecimal output but curiously enough not in the actual data (I write out the variables)
>
> Very odd. I must be looking at it in a wrong way.
>
>

Can you post both the Fortran program and the C program?

Louis

On Friday, September 2, 2022 at 7:38:46 AM UTC+2, Louis Krupp wrote:
> On 9/1/2022 11:24 PM, Arjen Markus wrote:
> > On Thursday, September 1, 2022 at 10:39:32 AM UTC+2, Arjen Markus wrote:
> >
> >> I tried to view the result of this Fortran program with ghostscript and got rubbish. I also tried the result of the C program and that showed the picture I expected. So there is something different, possibly caused by the formatting. But I am no expert in PostScript.
> >>
> > Since I get a completely contorted picture with the Fortran program I am trying to figure out what is going on. The results are identical between Intel Fortran oneAPI and gfortran, but differ greatly when using the C version.
> >
> > I have tried analysing the differences between the output from the Fortran program and the C program. So far I can find no cause, but:
> > - The output file from the C program is roughly 54000 bytes larger
> > - The number of lines is the same (I used a similar modification as John, but by filling a string of 60 characters)
> > - The lengths of the lines are equal
> > - The differences occur in the hexadecimal output but curiously enough not in the actual data (I write out the variables)
> >
> > Very odd. I must be looking at it in a wrong way.
> >
> >
> Can you post both the Fortran program and the C program?
>
> Louis

Here is the C code with the extra output:

#include <stdio.h>
#include <math.h>

#define MULT(x,y,a,b,t) ((t)=(x)*(a)-(y)*(b),(y)=(x)*(b)+(y)*(a),(x)=(t))
#define DIV(x,y,a,b,t) (MULT(x,y,a,-b,t),x /= SQ(a,b),y /= SQ(a,b))
#define SQ(x,y) ((x)*(x)+(y)*(y))
#define pi 3.141592654
#define R 150

#define CONST 1.905e-3
#define ROWS 900
#define COLS 900
#define R2 .16

main() {
FILE *out;
double x,y,zr,zi,wr,wi,dr,di,angle,t;
int i,j,gray,ix,iy;
static char dots[]={0x80,0x40,0x20,0x10,0x08,0x04,0x02,0x01};
static char *ps[]={"%%!\n",
"/str 256 string def\n",
"/go { image showpage } def\n",
"0 720 translate 72 72 scale\n",
"%d %d 8 [%d 0 0 -%d 0 %d]\n",
"{currentfile str readhexstring pop}\n",
"go\n"};

out = fopen("rootc.out", "w");

for(i=0;i<sizeof(ps)/sizeof(char*);i++) printf(ps[i],ROWS,COLS,R,R,R);

j=0;
for(ix=-ROWS/2;ix<ROWS-ROWS/2;ix++)
for(iy=-COLS/2;iy<COLS-COLS/2;iy++) {
zr=ix*CONST;
zi=iy*CONST;
gray=65535;
for(i=1;i<=224;i++) {
wr=zr,wi=zi;
MULT(wr,wi,zr,zi,t);
MULT(zr,zi,wr,wi,t);
zr *=2, zi *= 2;
zr += 1;
DIV(zr,zi,wr,wi,t);
zr /= 3; zi /= 3;
dr=zr-1,di=zi;
if(SQ(dr,di)<R2) break;
dr=zr+0.5,di=zi+0.866;
if(SQ(dr,di)<R2) break;
dr=zr+0.5,di=zi-0.866;
if(SQ(dr,di)<R2) break;
} if(i>224) angle=atan2(0.,-1.);
else angle=atan2(di,dr);
#if 0
fprintf(stderr,"%d %d %f %d\n",ix,iy,angle,i);
#endif
printf("%02x",(int)(fabs(angle)/.012272));
if(++j%30==0) putchar('\n');

fprintf(out, "%5d%5d %5d %16.8lf%16.8lf%16.8lf%16.8lf%16.8lf%16.8lf\n", ix, iy, (int)(fabs(angle)/.012272), ix*CONST, iy*CONST, zr, zi, dr, di);
} }

and here is the Fortran code:

PROGRAM ROOT3
REAL*8 PI/3.141592653589793d0/, CONST/1.905d-3/, R2R/.4d0/
INTEGER ROWS/900/, COLS/900/, R/150/
REAL*8 ANGLE
COMPLEX*16 Z, W, D
INTEGER I, J, GRAY, IX, IY
CHARACTER*50 PS(7)
character(len=60) :: output
integer :: pos
DATA PS/"%!",
* "/str 256 string def",
* "/go { image showpage } def",
* "0 720 translate 72 72 scale",
* "(I6, I6, ' 8 [', I6, ' 0 0 ', I6, ' 0 ', I6, ']')",
* "{currentfile str readhexstring pop}",
* "go"/
DO I=1,4
WRITE(*,'(A)') PS(I)
END DO
WRITE(*,PS(5)) ROWS, COLS, R, -R, R
DO I=6,7
WRITE(*,'(A)') PS(I)
END DO

pos = 1

DO IX=-ROWS/2, ROWS-ROWS/2-1
DO IY=-COLS/2, COLS-COLS/2-1
Z=CMPLX(IX*CONST, IY*CONST, KIND(1.D0))
DO I=1,224
W=Z*Z
Z=Z*W*2+1
Z=Z/W/3
D=Z-1
IF(ABS(D)<R2R) EXIT
D=Z+CMPLX(0.5D0, 0.866D0, KIND(1.D0))
IF(ABS(D)<R2R) EXIT
D=Z+CMPLX(0.5D0, -0.866D0, KIND(1.D0))
IF(ABS(D)<R2R) EXIT
END DO
IF(I>224) THEN
ANGLE=ATAN2(0.d0,1.d0)
ELSE
ANGLE=ATAN2(AIMAG(D), REAL(D))
ENDIF
! WRITE(*,'(Z2.2)') INT(ABS(ANGLE/0.012272))
WRITE(output(pos:pos+1),'(Z2.2)')
& INT(ABS(ANGLE/0.012272d0))
if ( pos < 60 ) then
pos = pos + 2
ELSE
pos = 1
write(*,'(a)') output
ENDIF
write(88,'(2i5,x,i5,x,10f16.8)') ix, iy,
& INT(ABS(ANGLE/0.012272d0)),
& ix*const, iy*const, z, d
END DO
END DO
END

The extra files are the same in size. by the way, minimal differences due to rounding, whereas the PS files are 1620166 bytes for the Fortran version and 1674155 for the C version.

Regards,

Arjen

On Friday, September 2, 2022 at 8:47:07 AM UTC+2, Arjen Markus wrote:
....

Unfortunately, all leading spaces are removed again :(

Regards,

Arjen

On 9/1/2022 4:10 PM, John wrote:
> Well, that worked for a long time, but google groups at least still shows no indenting. Oops. That thing goes out of it's way to mess up any kind of formatting, etc.

The indentation is fine when I view your post (with the Fortran code)
using Thunderbird. Copy-paste-compile works, as well.

Anyone using a news reader in combination with a Usenet service should
see the code fine, as well.

-- CyrMag

Test for pos < 59 instead of 60 because incrementing by two and I think the Fortran code would work. The example is divisible by 30 so the output line always gets filled, but if the size were changed to something not divisible by 30 the output would be full of the previous pass and not get written; but is OK as it is with the current size.

On Friday, September 2, 2022 at 4:04:11 PM UTC+2, John wrote:
> Test for pos < 59 instead of 60 because incrementing by two and I think the Fortran code would work. The example is divisible by 30 so the output line always gets filled, but if the size were changed to something not divisible by 30 the output would be full of the previous pass and not get written; but is OK as it is with the current size.

Indeed, with the revised test I get the expected intricate picture! Thanks. There is a single pixel that has a different value in the Fortran version than in the C version, but I doubt anyone would see that.

Regards,

Arjen