C Program hexpet16.f to create any number of concentric circles C of distorted hexagons, each with corners on its own circum-circle, C with a circular boundary around more rings of hexagons. C The centres of the hexagons in the outer circle lie C at a radius R = 1.0 from the centre of the pattern. C The radius of the circum-circle of the outer hexagons is HR C The scaling factor from one circle to the next is Q C The position angles of the corners Th(1) to Th(6) are measured C clockwise from the radial line from the centre of the pattern C through the centre of the circum-circle. Position angles of C the centres of the circum-circles in the second ring are Th(7) C and Th(8). C Implicit Real*8(A-H,O-Z) Real*4 ymin, ymax, zmin, zmax, ycent, zcent, rad, Y, Z, W, X Dimension U(20), V(20), Y(1999), Z(1999), Th(10), R(20), 1 S(20), W(6), X(6) character infile*80, outfile*80, message*80 Zero = 0.0D0 One = 1.0D0 Pi = 3.141592653589793D0 Radian = Pi/1.8D2 ATPi = Pi/3.0D0 TTPi = 2.0D0*Pi/3.0D0 200 write(6,*) ' Name the input file . . .dat' read(5,'(a)') infile open(7,file=infile,status='old',err=200) 210 write(6,*) ' and the output file . . .out' read(5,'(a)') outfile open(8,file=outfile,status='new',err=210) 118 Format(4D10.5,2I5) Read (7,120) Step, Nring 120 Format(D7.1,I5) C HR and Q are fine-tuned to allow the iterative process to converge HR = 3.62D0/Step - 5.12D0/(Step**3) HR = (Dsin(Pi/Step)/Dsin(Pi/3.0D0)) A = Dsin((Pi/6.0D0) - (Pi/(2.0D0*Step))) B = Dsin((Pi/6.0D0) + (Pi/(2.0D0*Step))) Q = A/B C Suppress diagnostic print-out: C Write (8,68) Q, HR 68 Format(1H ,'Trial values of Q and HR',2F20.15) Trad = -0.1D0 Th(1) = 0.0D0 Th(2) = ATPi + Pi/Step Th(3) = ATPi + Th(2) Th(4) = Pi Th(6) = Pi + Pi - ATPi - Pi/Step Th(5) = Th(6) - ATPi DTh = 0.01D0 Icount = 0 Iflag = 0 Jflag = 0 100 Continue If(Iflag.NE.0) Go To 1 C Read frame size of new pgplot window. Read (7,118) ymin, ymax, zmin, zmax Call pgbegin (0,'?',1,1) ymin = Sngl(ymin) ymax = Sngl(ymax) zmin = Sngl(zmin) zmax = Sngl(zmax) IP = Ifix(Sngl(Step)) Call pgenv (ymin, ymax, zmin, zmax, 0, -2) C Write (message,222) Step 222 Format(1H ,'Concentric rings of',F6.1,' hexagons') C Call pglabel (' ',' ',message) Call pgsfs (2) Call pgsls (1) Call pgslw (2) 1 Th(9) = 0.0D0 Th(8) = Pi/Step 3 Th(7) = -Pi/Step C Th(9) is the position angle of a hexagon in the outer circle. C Th(7) and Th(8) are position angles of the centres of hexagons C adjacent to the first one. Q is the factor by which hexagons C shrink on going from one circle to the next one inside C of the start of the next spiral. CR is the distance of the C centre of the first hexagon from the centre of the pattern. C CR1 is the distance of the second and third hexagons C from the centre of the pattern. HR1 = HR*Q CR = Trad If(Trad.GT.Zero) Go To 14 CR = One 14 CR1 = CR*Q Do 4 I = 1,6 U(I) = CR*Dsin(Th(9)) + HR*Dsin(Th(I)) V(I) = CR*Dcos(Th(9)) + HR*Dcos(Th(I)) U(I+6) = CR1*Dsin(Th(7)) + HR1*Dsin(Th(7)+Th(I)) V(I+6) = CR1*Dcos(Th(7)) + HR1*Dcos(Th(7)+Th(I)) U(I+12) = CR1*Dsin(Th(8)) + HR1*Dsin(Th(8)+Th(I)) V(I+12) = CR1*Dcos(Th(8)) + HR1*Dcos(Th(8)+Th(I)) 4 Continue Do 102 I = 1,18 Y(I) = Sngl(U(I)) Z(I) = Sngl(V(I)) 102 Continue If(Iflag.NE.0) Go To 103 Do 104 I = 1,6 Y(I) = Sngl(U(I)) Z(I) = Sngl(V(I)) 104 Continue Y(7) = Sngl(U(1)) Z(7) = Sngl(V(1)) Do 105 I = 1,6 Y(I) = Sngl(U(I+6)) Z(I) = Sngl(V(I+6)) 105 Continue Y(7) = Sngl(U(7)) Z(7) = Sngl(V(7)) Do 106 I = 1,6 Y(I) = Sngl(U(I+12)) Z(I) = Sngl(V(I+12)) 106 Continue Y(7) = Sngl(U(13)) Z(7) = Sngl(V(13)) Y(1) = Sngl(CR*Sin(Th(9))) Z(1) = Sngl(CR*Cos(Th(9))) Y(2) = Sngl(CR1*Sin(Th(7))) Z(2) = Sngl(CR1*Cos(Th(7))) Y(3) = Sngl(CR1*Sin(Th(8))) Z(3) = Sngl(CR1*Cos(Th(8))) Ycent = Sngl(Y(1)) Zcent = Sngl(Z(1)) rad = Sngl(HR) Ycent = Sngl(Y(2)) Zcent = Sngl(Z(2)) rad = Sngl(HR1) Ycent = Sngl(Y(3)) Zcent = Sngl(Z(3)) rad = Sngl(HR1) 103 Continue C Diagnostic printout: C Write(8,101) 101 Format(1H ) C Write(8,5) Q, Step 5 Format(1H ,'Q, Step ',2F20.15) C Write(8,6) CR, CR1 6 Format(1H ,'CR, CR1 ',2F20.15) C Write(8,7) HR, HR1 7 Format(1H ,'HR, HR1 ',2F20.15) C Write(8,101) Do 9 I = 1,18 R(I) = Dsqrt(U(I)*U(I) + V(I)*V(I)) C Write(8,8) I, U(I), V(I), R(I) 8 Format(1H ,I5,3F15.9) 9 Continue Do 10 I = 1,5 S(I) = Dsqrt((U(I)-U(I+1))**2 + (V(I)-V(I+1))**2) S(I+6) = Dsqrt((U(I+6)-U(I+7))**2 + (V(I+6)-V(I+7))**2) S(I+12) = Dsqrt((U(I+12)-U(I+13))**2 + (V(I+12)-V(I+13))**2) 10 Continue S(6) = Dsqrt((U(6)-U(1))**2 + (V(6)-V(1))**2) S(12) = Dsqrt((U(12)-U(7))**2 + (V(12)-V(7))**2) S(18) = Dsqrt((U(18)-U(13))**2 + (V(18)-V(13))**2) C Write(8,101) C Write(8,11) 11 Format(1H ,' Sides clockwise from corner no.') C Write(8,101) Do 13 I = 1,18 C Write(8,12) I, S(I) 12 Format(1H ,I5,F15.9) 13 Continue 23 If(V(4).GT.V(8)) Go To 24 Q = Q - 0.1D0*DTh Go To 25 24 Q = Q + 0.1D0*DTh 25 If((U(17)+U(18)).GT.(U(8)+U(9))) Go To 26 HR = HR - 0.1D0*DTh Go To 27 26 HR = HR + 0.1D0*DTh 27 Continue If (Iflag.NE.0) Go To 36 36 Iflag = Iflag + 1 Icount = Icount + 1 If(Iflag.LT.50) Go To 37 C Iflag = 0 37 DTh = DTh*0.90D0 C Suppress diagnostic print-out: C Write (8,47) Iflag, Icount, Q, HR 47 Format(1H ,2I5,2F20.15) If(DTh.GT.3.0D-16) Go To 100 Write (8,101) Write (8,5) Q, Step Write (8,101) Write (8,6) CR, CR1 Write (8,101) Write (8,7) HR, HR1 Write (8,101) Write (8,96) 96 Format(1H ,' I U V R 1 S') Write (8,101) Do 99 I = 1,18 Write (8,97) I, U(I), V(I), R(I), S(I) 97 Format(1H ,I4,4F17.13) 99 Continue Write (8,101) Write (8,95) (Th(I), I = 1,6) Temp1 = Th(2)/Radian Temp2 = Th(3)/Radian Temp3 = Th(5)/Radian Temp4 = Th(6)/Radian Write (8,98) Temp1, Temp2, Temp3, Temp4 98 Format(1H ,'P.A. of corners ',4F14.9) Temp1 = Th(6) - Th(5) + Th(3) - Th(2) Temp2 = Th(5) - Th(4) + Th(2) - Th(1) Temp3 = Th(4) - Th(3) + Th(1) - Th(6) + Pi + Pi Write (8,95) Temp1, Temp2, Temp3 95 Format(1H ,3F20.16) Write (8,101) Write (8,70) 70 Format(1H ,'Entries in the third column should all be zeros.') Temp = R(6) Temp1 = Temp - R(2) Write (8,94) Temp, R(2), Temp1 94 Format(1H ,'R(6) and R(2) ',3F17.13) Temp = R(5) Temp1 = Temp - R(3) Write (8,93) Temp, R(3), Temp1 93 Format(1H ,'R(5) and R(3) ',3F17.13) Temp = R(7) Temp1 = Temp - R(13) Write (8,92) Temp, R(13), Temp1 92 Format(1H ,'R(7) and R(13) ',3F17.13) Temp = R(12) Temp1 = Temp - R(8) Write (8,91) Temp, R(8), Temp1 91 Format(1H ,'R(12) and R(8) ',3F17.13) Temp = R(14) Temp1 = Temp - R(18) Write (8,90) Temp, R(18), Temp1 90 Format(1H ,'R(14) and R(18) ',3F17.13) Temp1 = R(4) - R(8) Write (8,89) R(4), R(8), Temp1 89 Format(1H ,'R(4) and R(8) ',3F17.13) Temp = R(11) Temp1 = Temp - R(17) Write (8,88) Temp, R(17), Temp1 88 Format(1H ,'R(11) and R(17) ',3F17.13) Temp = R(15) Temp1 = Temp - R(9) Write (8,87) Temp, R(9), Temp1 87 Format(1H ,'R(15) and R(9) ',3F17.13) Temp = R(10) Temp1 = Temp - R(16) Write (8,86) Temp, R(16), Temp1 86 Format(1H ,'R(10) and R(16) ',3F17.13) Write (8,101) Write (8,71) 71 Format(1H ,'Entries in all three columns should all be zeros.') Temp1 = U(5)-U(7) Temp2 = V(5)-V(7) Temp3 = Sqrt(Temp1*Temp1 + Temp2*Temp2) Write (8,85) Temp1, Temp2, Temp3 85 Format(1H ,'dU, dV, Total at 5/7',3F17.13) Temp1 = U(3)-U(13) Temp2 = V(3)-V(13) Temp3 = Sqrt(Temp1*Temp1 + Temp2*Temp2) Write (8,84) Temp1, Temp2, Temp3 84 Format(1H ,'dU, dV, Total at 3/13',3F17.13) Temp1 = U(4)-U(8) Temp2 = V(4)-V(8) Temp3 = Sqrt(Temp1*Temp1 + Temp2*Temp2) Write (8,83) Temp1, Temp2, Temp3 83 Format(1H ,'dU, dV, Total at 4/8',3F17.13) Temp1 = U(8)-U(18) Temp2 = V(8)-V(18) Temp3 = Sqrt(Temp1*Temp1 + Temp2*Temp2) Write (8,82) Temp1, Temp2, Temp3 82 Format(1H ,'dU, dV, Total at 8/18',3F17.13) Temp1 = U(9)-U(17) Temp2 = V(9)-V(17) Temp3 = Sqrt(Temp1*Temp1 + Temp2*Temp2) Write (8,81) Temp1, Temp2, Temp3 81 Format(1H ,'dU, dV, Total at 9/17',3F17.13) Write (8,101) Write (8,60) 60 Format(1H ,'Entries in the third column should all be zeros.') Temp = S(6) Temp1 = Temp - S(1) Write (8,80) Temp, S(1), Temp1 80 Format(1H ,'S(6) and S(1) ',3F17.13) Temp = S(5) Temp1 = Temp - S(2) Write (8,79) Temp, S(2), Temp1 79 Format(1H ,'S(5) and S(2) ',3F17.13) Temp = S(4) Temp1 = Temp - S(3) Write (8,78) Temp, S(3), Temp1 78 Format(1H ,'S(4) and S(3) ',3F17.13) Temp = S(1)*Q Temp1 = Temp - S(4) Write (8,77) Temp, S(4), Temp1 77 Format(1H ,'S(1)*Q and S(4) ',3F17.13) Temp = S(6)*Q Temp1 = Temp - S(3) Write (8,76) Temp, S(3), Temp1 76 Format(1H ,'S(6)*Q and S(3) ',3F17.13) Temp1 = S(4) - S(7) Write (8,75) S(4), S(7), Temp1 75 Format(1H ,'S(4) and S(7) ',3F17.13) Temp1 = S(3) - S(18) Write (8,74) S(3), S(18), Temp1 74 Format(1H ,'S(3) and S(18) ',3F17.13) Temp1 = S(8) - S(17) Write (8,73) S(8), S(17), Temp1 73 Format(1H ,'S(8) and S(17) ',3F17.13) C Several lines removed from program here, but it still works. Angle1 = -Pi/Step Angle2 = 0.0D0 Write (8,101) Write (8,72) Angle1, Angle2 72 Format(1H ,'Starting angles for circular pattern are',2F17.13) Dangle = 2.0D0*Pi/Step Diag1 = R(1) - R(4) Diag2 = R(2)*R(2) + R(5)*R(5) - 2.0D0*R(2)*R(5)*Dcos(Dangle) Diag2 = Dsqrt(Diag2) Write (8,59) Diag1, Diag2 59 Format(1H ,'Radial and oblique diags in outer ring ',2F17.13) IPP = IP + IP + 2 C Generate the outer ring of corners of hexagons Do 61 I = 1, IPP, 2 Y(I) = Sngl(R(6)*Dsin(Angle1)) Z(I) = Sngl(R(6)*Dcos(Angle1)) Y(I+1) = Sngl(R(1)*Dsin(Angle2)) Z(I+1) = Sngl(R(1)*Dcos(Angle2)) Angle1 = Angle1 + Dangle Angle2 = Angle2 + Dangle 61 Continue C Draw the outer ring of corners if Trad is -ve. If(Trad.GT.Zero) Go To 46 Call pgline (IPP,y,z) C Calculate the distance between corners 2 and 6, and print it. Write (8,101) Size = 2.0D0*R(6)*Dsin(Angle1) Write (8,43) Step, Size 43 Format(1H ,F6.1,' pieces, max. circumf. is ',F12.8) Go To 45 C If Trad is positive, reset the count of rings, and scale the C rings to fit 46 Lring = 0 Do 44 I = 1,IPP Y(I) = Sngl(Y(I)*Trad/R(6)) Z(I) = Sngl(Z(I)*Trad/R(6)) 44 Continue R(5) = R(5)*Trad/R(6) R(4) = R(4)*Trad/R(6) R(3) = R(3)*Trad/R(6) R(2) = R(2)*Trad/R(6) R(1) = R(1)*Trad/R(6) R(6) = Trad Angle1 = -Pi/Step Write (8,101) Size = - 2.0D0*R(6)*Dsin(Angle1) Write (8,43) Step, Size 45 Angle1 = -Pi/Step Angle2 = 0.0D0 C Generate second, fourth, etc. ring of corners Do 63 I = 1,IPP,2 Y(I+IPP) = Sngl(R(5)*Dsin(Angle1)) Z(I+IPP) = Sngl(R(5)*Dcos(Angle1)) Y(I+IPP+1) = Sngl(R(4)*Dsin(Angle2)) Z(I+IPP+1) = Sngl(R(4)*Dcos(Angle2)) Angle1 = Angle1 + Dangle Angle2 = Angle2 + Dangle 63 Continue C Draw radial lines, data transferred to arrays w and x. K = 1 W(1) = Sngl(Y(K)) X(1) = Sngl(Z(K)) W(2) = Sngl(Y(K+IPP)) X(2) = Sngl(Z(K+IPP)) Rrad = R(5) rad = Sngl(Rrad) Write (6,42) rad, R(4), R(5), Lring 42 Format(1H ,'Line 372 rad, R(4), R(5) =',3F9.5,I5) Call pgline (2,w,x) Do 62 K = 3,IPP,2 C Outline hexagons, corners in arrays w and x W(6) = Sngl(Y(K-1)) X(6) = Sngl(Z(K-1)) W(5) = Sngl(W(1)) X(5) = Sngl(X(1)) W(4) = Sngl(W(2)) X(4) = Sngl(X(2)) W(3) = Sngl(Y(K+IPP-1)) X(3) = Sngl(Z(K+IPP-1)) W(2) = Sngl(Y(K+IPP)) X(2) = Sngl(Z(K+IPP)) W(1) = Sngl(Y(K)) X(1) = Sngl(Z(K)) C Draw radial lines Call pgline (2,w,x) 62 Continue Lring = Lring + 1 If(Lring.GE.Nring) Go To 57 C Shift the second set of IPP points Do 56 K = 1,IPP Y(K) = Sngl(Y(K+IPP)) Z(K) = Sngl(Z(K+IPP)) 56 Continue Rzig = R(5) rad = Sngl(Rzig) Write (6,40) rad, R(4), R(5), Lring 40 Format(1H ,'Line 405 rad, R(4), R(5) =',3F9.5,I5) Jflag = Jflag + 2 55 Call pgline (IPP,y,z) 52 Angle1 = -Pi/Step Angle2 = 0.0D0 C Generate the third, fifth, etc. ring of corners R(1) = R(1)*Q*Q R(2) = R(2)*Q*Q R(4) = R(4)*Q*Q R(5) = R(5)*Q*Q Do 64 I = 1,IPP,2 Y(I+IPP) = Sngl(R(2)*Dsin(Angle1)) Z(I+IPP) = Sngl(R(2)*Dcos(Angle1)) Y(I+IPP+1) = Sngl(R(1)*Dsin(Angle2)) Z(I+IPP+1) = Sngl(R(1)*Dcos(Angle2)) Angle1 = Angle1 + Dangle Angle2 = Angle2 + Dangle 64 Continue C Draw radial lines, with data transferred to arrays w and x. K = 2 W(1) = Sngl(Y(K)) X(1) = Sngl(Z(K)) W(2) = Sngl(Y(K+IPP)) X(2) = Sngl(Z(K+IPP)) Rrad = R(5)/Q rad = Rrad Write (6,41) rad, R(4), R(5), Lring 41 Format(1H ,'Line 435 rad, R(4), R(5) =',3F9.5,I5) Jflag = Jflag + 2 53 Call pgline (2,w,x) Do 65 K = 4,IPP,2 C Outline hexagons, corners in arrays w and x W(6) = Sngl(Y(K-1)) X(6) = Sngl(Z(K-1)) W(5) = Sngl(W(1)) X(5) = Sngl(X(1)) W(4) = Sngl(W(2)) X(4) = Sngl(X(2)) W(3) = Sngl(Y(K+IPP-1)) X(3) = Sngl(Z(K+IPP-1)) W(2) = Sngl(Y(K+IPP)) X(2) = Sngl(Z(K+IPP)) W(1) = Sngl(Y(K)) X(1) = Sngl(Z(K)) C Draw radial lines Call pgline (2,w,x) 67 Continue 65 Continue Lring = Lring + 1 If(Lring.GE.Nring) Go To 57 50 Continue Do 66 K = 1,IPP Y(K) = Sngl(Y(K+IPP)) Z(K) = Sngl(Z(K+IPP)) 66 Continue Rzag = Sqrt(Y(2)*Y(2) + Z(2)*Z(2)) Call pgline(IPP,y,z) 54 Continue Go To 45 57 Continue Call pgcirc (0.0, 0.0, rad) Read (7,120) Step, Nring Write (6,69) rad, Trad, Rzag, R(4), R(5) 69 Format(1H ,'r, Tr, Rz, R(4,5) = ',5F9.5) Trad = R(5) C HR and Q are fine-tuned to allow the iterative process to converge HR = Trad*(3.62D0/Step - 5.12D0/(Step**3)) HR = Trad*((Dsin(Pi/Step))/(Dsin(Pi/3.0D0))) A = Dsin((Pi/6.0D0) - (Pi/(2.0D0*Step))) B = Dsin((Pi/6.0D0) + (Pi/(2.0D0*Step))) Q = A/B C Suppress diagnostic print-out: C Write (8,68) Q, HR Th(1) = Zero Th(2) = ATPi + Pi/Step Th(3) = ATPi + Th(2) Th(4) = Pi Th(6) = Pi + Pi - ATPi - Pi/Step Th(5) = Th(6) - ATPi DTh = 0.01D0 Icount = 0 Iflag = 0 Jflag = 0 C XX Write (8,48) Q, HR 48 Format(1H ,'New cycle',2F20.15) If (Step.GT.Zero) Go To 1 Write (8,51) rad 51 Format(1H ,'Radius of central monolith is ',F9.6) rad = Sngl(0.35*rad) Call pgcirc (0.0, 0.0, rad) 58 Call pgend STOP End