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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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');
}
}

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

Unfortunately, all leading spaces are removed again :(

Regards,

Arjen

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)

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

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)