PROGRAM hyper IMPLICIT NONE REAL(KIND=4), DIMENSION(:,:,:), ALLOCATABLE :: realdata REAL(KIND=4) :: x, y, z, dx, minx, maxx, arg INTEGER(KIND=2), DIMENSION(:,:,:), ALLOCATABLE :: odata INTEGER(KIND=2) :: NX, NY, NZ INTEGER(KIND=4) :: J, K, L, slice CHARACTER(LEN=60) :: ofile !Set input parameters (size of grid) !----------------------------------------------------------------------------------! minx = -2.0 maxx = 2.0 dx = 0.02 !----------------------------------------------------------------------------------! !Determine size of grid - crude! x = minx NX = 0 DO WHILE (x <= maxx) x = x + dx NX = NX + 1 ENDDO NY = NX NZ = NX write(*,*) "Grid Spacing: ", NX, NY, NZ !Allocate data ALLOCATE(realdata(NX,NY,NZ)) ALLOCATE(odata(NX,NY,NZ)) !Generate the volume data x = minx J = 1 DO WHILE (J <= NX) y = minx K = 1 DO WHILE (K <= NY) z = minx L = 1 DO WHILE (L <= NZ) !sphere !----------------------------------------------------------------------------! ofile = 'sphere.df3' realdata(J,K,L) = x**2 + y**2 + z**2 !----------------------------------------------------------------------------! !5-hole bretzel !----------------------------------------------------------------------------! !ofile = 'bretzel.df3' !arg = (x*x+y*y/4. - 1.)*(x*x/4.0+y*y-1.0) !realdata(J,K,L) = arg**2 + z*z -0.1 !----------------------------------------------------------------------------! !Goursat ! works best with dx = 0.02, minx = -2.0, maxx = 2.0 !----------------------------------------------------------------------------! !ofile = 'goursat.df3' !realdata(J,K,L) = x**4 + y**4 + z**4 + (-0.27)*((x**2 + y**2 + z**2)**2) & !& + (-0.5)*(x**2 + y**2 + z**2) -1.0 !----------------------------------------------------------------------------! !Cassini ! works with dx = 0.02, minx = -2.0, maxx = 2.0 !----------------------------------------------------------------------------! !ofile = 'cassini.df3' !arg = 0.3 !realdata(J,K,L) = ((x-arg)**2 + z**2)*((x+arg)**2 + z**2) - y**4 !Barth ! works with dx = 0.02, minx = -2.0, maxx = 2.0 !----------------------------------------------------------------------------! !ofile = 'barth.df3' !arg = 1.6180339887 !realdata(J,K,L) = 4.0*((arg**2)*(x**2)-y**2)*((arg**2)*(y**2)-z**2)* & !& ((arg**2)*(z**2)-x**2) - (1. + 2.*arg)* & !& ((x**2 + y**2 + z**2 - (0.9)**2)**2)*(0.9)**2 z = z + dx L = L + 1 ENDDO y = y + dx K = K + 1 ENDDO x = x + dx J = J + 1 ENDDO !Let's write out integer values, scaled to 2-bytes !First scale real values to be greater than zero IF (MINVAL(realdata(:,:,:)) < 0) THEN realdata(:,:,:) = realdata(:,:,:) + ABS(MINVAL(realdata(:,:,:))) ENDIF !Now maximize the dynamic range of the values to the maximum 2-byte integer size realdata(:,:,:) = (realdata(:,:,:)/MAXVAL(realdata(:,:,:)))*65534 !Now put these values into an array of 2-byte integers odata(:,:,:) = 0 odata(:,:,:) = NINT(realdata(:,:,:)) !Now check to make sure we are still within the allowable range of ! 2-byte integers write(*,*) "Min/Max Integers: ", MINVAL(odata(:,:,:)), MAXVAL(odata(:,:,:)) IF (MINVAL(odata(:,:,:)) < -32767 .OR. MAXVAL(odata(:,:,:)) > 32767) THEN write(*,*) "WARNING: Check your integer range!" ENDIF !Write out the density file! CALL writedensfile(NX,NY,NZ,odata,ofile) !Write out a test slice on the interval 0 to 1 and scaled from 0 to 1 !realdata(:,:,:) = realdata(:,:,:)/MAXVAL(realdata(:,:,:)) !OPEN(UNIT=1,FILE='testslice.xy') !x = 0.0 !slice = (NX-1)/2 + 1 !DO J=1,NX ! y = 0.0 ! DO K=1,NY ! write(1,*) x, y, realdata(J,K,slice) ! y = y + dx ! ENDDO !x = x + dx !ENDDO !CLOSE(1) END PROGRAM hyper !--------------------------------------------------------------------------------------------------------! ! WRITEDENSFILE ! ! Subroutine to write out .df3 format file (density file) for use with POVRay ! Note, this subroutine is designed to write out 2-byte integers for the density values. ! !--------------------------------------------------------------------------------------------------------! SUBROUTINE writedensfile(NX,NY,NZ,outdata,ofile) IMPLICIT NONE INTEGER(KIND=2), INTENT(IN) :: NX, NY, NZ !array dimensions INTEGER(KIND=2), DIMENSION(NX,NY,NZ), INTENT(IN) :: outdata !density values to write out INTEGER(KIND=4) :: x, y, z !loop variables INTEGER(KIND=4) :: counter !record counter for direct access output INTEGER(KIND=2) :: bitx, bity, bitz, bitd !test bit size CHARACTER(LEN=60), INTENT(IN) :: ofile !name of output file to write to ! First test bit size, header MUST always be written as 2-byte integers ! In this subroutine I use KIND=2 to indicate 2-byte integers. But, ! note that some compilers use a sequential KIND ordering, and KIND=2 ! may not always indicate 2-byte integers. This is intended to check ! that possibility. bitx = bit_size(NX) bity = bit_size(NY) bitz = bit_size(NZ) bitd = bit_size(outdata) IF (bitx == 16 .AND. bity == 16 .AND. bitz == 16 .AND. bitd == 16) THEN write(*,*) "Writing density file: '", TRIM(ADJUSTL(ofile)), "'..." write(*,*) "Expected file size is: ", 2*(NX*NY*NZ) + 6, " (bytes)" ELSE write(*,*) "ERROR: INTEGER bit size does not equal 16. Exiting Now..." write(*,*) "Bit Size of NX: ", bitx write(*,*) "Bit Size of NY: ", bity write(*,*) "Bit Size of NZ: ", bitz write(*,*) "Bit Size of output data: ", bitd STOP ENDIF ! First write out 6-byte header with array dimensions ! Output file is binary and MUST be BIG_ENDIAN OPEN(UNIT=1,FILE=ofile,ACCESS='DIRECT',FORM='UNFORMATTED',RECL=6,STATUS='REPLACE',CONVERT='BIG_ENDIAN') WRITE(1,rec=1) NX, NY, NZ CLOSE(1) ! Now write out density values counter = 4 OPEN(UNIT=1,FILE=ofile,ACCESS='DIRECT',FORM='UNFORMATTED',RECL=2,CONVERT='BIG_ENDIAN') DO z=1,NZ DO y=1,NY DO x=1,NX write(1,rec=counter) outdata(x,y,z) counter = counter + 1 ENDDO ENDDO ENDDO CLOSE(1) END SUBROUTINE writedensfile !--------------------------------------------------------------------------------------------------------!