! File: fire2.f90
! Public domain 2006 James Van Buskirk
module bitmap_types
   implicit none
   integer, parameter :: ik1 = selected_int_kind(2)
   integer, parameter :: ik2 = selected_int_kind(4)
   integer, parameter :: ik4 = selected_int_kind(9)
   integer, parameter :: ik8 = selected_int_kind(18)
   integer, parameter :: sp = selected_real_kind(6,30)
   integer, parameter :: dp = selected_real_kind(15,300)

   type BITMAPFILEHEADER
      SEQUENCE
! df & ifort pack SEQUENCE types by default; lf95 & g95 do not
!      integer(ik2) bfType       ! signature word "BM" or Z'4d42'
!      integer(ik4) bfSize       ! entire size of file
!      integer(ik2) bfReserved1  ! must be zero
!      integer(ik2) bfReserved2  ! must be zero
!      integer(ik4) bfOffsetBits ! offset in file of DIB pixel bits
      integer(ik2) bfType          ! signature word "BM" or Z'4d42'
      integer(ik2) bfSize(2)       ! entire size of file
      integer(ik2) bfReserved1     ! must be zero
      integer(ik2) bfReserved2     ! must be zero
      integer(ik2) bfOffsetBits(2) ! offset in file of DIB pixel bits
   end type BITMAPFILEHEADER

   type RGBQUAD
      SEQUENCE
      integer(ik1) rgbBlue     ! blue level
      integer(ik1) rgbGreen    ! green level
      integer(ik1) rgbRed      ! red level
      integer(ik1) rgbReserved ! = 0
   end type RGBQUAD

   type CIEXYZ
      SEQUENCE
      integer(ik4) ciexyzX
      integer(ik4) ciexyzY
      integer(ik4) ciexyzZ
   end type CIEXYZ

   type CIEXYZTRIPLE
      SEQUENCE
      type(CIEXYZ) ciexyzRed
      type(CIEXYZ) ciexyzGreen
      type(CIEXYZ) ciexyzBlue
   end type CIEXYZTRIPLE

   type BITMAPV5HEADER
      SEQUENCE
      integer(ik4) bV5Size            ! size of the structure = 120 (124?)
      integer(ik4) bV5Width           ! width of the image in pixels
      integer(ik4) bV5Height          ! height of the image in pixels
      integer(ik2) bV5Planes          ! = 1
      integer(ik2) bV5BitCount        ! bits per pixel (1, 4, 8, 16, 24, or 32)
      integer(ik4) bV5Compression     ! compression code
      integer(ik4) bV5SizeImage       ! number of bytes in image
      integer(ik4) bV5XPelsPerMeter   ! horizontal resolution
      integer(ik4) bV5YPelsPerMeter   ! vertical resolution
      integer(ik4) bV5ClrUsed         ! number of colors used
      integer(ik4) bV5ClrImportant    ! number of important colors
      integer(ik4) bV5RedMask         ! Red color mask
      integer(ik4) bV5GreenMask       ! Green color mask
      integer(ik4) bV5BlueMask        ! Blue color mask
      integer(ik4) bV5AlphaMask       ! Alpha mask
      integer(ik4) bV5CSType          ! color space type
      type(CIEXYZTRIPLE) bV5EndPoints ! XYZ values
      integer(ik4) bV5GammaRed        ! Red gamma value
      integer(ik4) bV5GammaGreen      ! Green gamma value
      integer(ik4) bV5GammaBlue       ! Blue gamma value
      integer(ik4) bV5Intent          ! rendering intent
      integer(ik4) bV5ProfileData     ! profile data or filename
      integer(ik4) bV5ProfileSize     ! size of embedded data or filename
      integer(ik4) bV5Reserved
   end type BITMAPV5HEADER
end module bitmap_types

module chk_ftn95
   implicit none
   private
   public is_ftn95
   interface operator(.NEQV.)
      module procedure my_neqv
   end interface operator(.NEQV.)
   contains
      elemental function my_neqv(a,b)
         logical my_neqv
         integer, intent(in) :: a, b

         my_neqv = (a == 0 .AND. b /= 0) .OR. &
                   (a /= 0 .AND. b == 0)
      end function my_neqv

      function is_ftn95()
         logical is_ftn95
         integer a, b

         a = 2
         b = 4
         if(a .NEQV. b) then
            is_ftn95 = .TRUE.
         else
            is_ftn95 = .FALSE.
         end if
      end function is_ftn95
end module chk_ftn95

module chk_g95
   implicit none
   private
   public is_g95
   type t
      integer x
   end type t
   interface assignment(=)
      module procedure assign_t
   end interface assignment(=)
   contains
      function is_g95()
         logical is_g95
         type(t) y(3)

         y%x = (/1,2,3/)
         y = y((/2,3,1/))
         is_g95 = y(3)%x == 1
      end function is_g95

      elemental subroutine assign_t(lhs,rhs)
         type(t), intent(in) :: rhs
         type(t), intent(out) :: lhs

         lhs%x = rhs%x
      end subroutine assign_t
end module chk_g95

program fire
   use bitmap_types
   use chk_ftn95
   use chk_g95
   implicit none
   type(BITMAPFILEHEADER) bmp_file_header
   type(RGBQUAD), allocatable :: palette(:)
   type(BITMAPV5HEADER) bmp_header
   integer RowLength, ColHeight
   integer(ik1), allocatable :: image(:,:)
   integer bfOffsetBits
   integer bfSize
   integer i, j
   integer PelsPerMeter
   real(dp) lambda
   integer x0, y0
   real(dp) f
   real(dp) r
   real(dp) width
   real(dp) height
   character(20) access, form
   integer DPI
   character(80) flt_err
   character(80) line
   integer line_size
   logical dither, apodize
   real(dp) harvest(2)
   real(dp) r0
   real(dp) sigma
   real(dp) r12
   real(dp) r1

! *************************************
! Create *.BMP file header
! File type code
   bmp_file_header%bfType = transfer('BM',bmp_file_header%bfType)
! (Done later) store entire size of file
   bmp_file_header%bfSize = transfer(0,bmp_file_header%bfSize)
! Reserved words alwys 0
   bmp_file_header%bfReserved1 = 0_ik2
   bmp_file_header%bfReserved2 = 0_ik2
! (Done later) store offset in file of DIB pixel bits
   bmp_file_header%bfOffsetBits = transfer(0,bmp_file_header%bfOffsetBits)

! *************************************
! Create *.BMP info header
   bmp_header%bV5Size = size(transfer(bmp_header,(/1_ik1/)))
   write(*,'(a)',advance='no') 'Enter resolution in dots/inch:> '
   read(*,*) DPI
   if(DPI <= 0) then
      write(*,'(a,i0,a)') 'Invalid resolution: ',DPI,' DPI'
      write(*,'(a)') 'Program aborting'
      stop
   end if
   write(*,'(a)',advance='no') 'Enter image width in inches:> '
   read(*,*) width
   if(width <= 0) then
      write(flt_err,*) width
      write(*,'(3a)') 'Invalid width: ',trim(adjustl(flt_err)),' inches'
      write(*,'(a)') 'Program aborting'
      stop
   end if
   write(*,'(a)',advance='no') 'Enter image height in inches:> '
   read(*,*) height
   if(height <= 0) then
      write(flt_err,*) height
      write(*,'(3a)') 'Invalid height: ',trim(adjustl(flt_err)),' inches'
      write(*,'(a)') 'Program aborting'
      stop
   end if
! (Done later) set image width in pixels
   bmp_header%bV5Width = 0
! (Done later) set image height in pixels
   bmp_header%bV5Height = 0
   bmp_header%bV5Planes = 1
! 1 bit per pixel
   bmp_header%bV5BitCount = 1
! Compression code = 0: no compression
   bmp_header%bV5Compression = 0
! (Done later) set image size in bytes
   bmp_header%bV5SizeImage = 0
! Horizontal resolution obtained above
   PelsPerMeter = nint(DPI/0.0254_dp)
   bmp_header%bV5XPelsPerMeter = PelsPerMeter
! Vertical resolution also obtained above
   bmp_header%bV5YPelsPerMeter = PelsPerMeter
! Use full palette
   if(bmp_header%bV5BitCount < 16) then
      bmp_header%bV5ClrUsed = 2**bmp_header%bV5BitCount
      bmp_header%bV5ClrImportant = bmp_header%bV5ClrUsed
   else
      bmp_header%bV5ClrUsed = 0
      bmp_header%bV5ClrImportant = 0
   end if
!! All colors important
!   bmp_header%bV5ClrImportant = 0
! We won't worry about the next 4 fields for now
   bmp_header%bV5RedMask   = 0
   bmp_header%bV5GreenMask = 0
   bmp_header%bV5BlueMask  = 0
   bmp_header%bV5AlphaMask = 0
   bmp_header%bV5CSType = 2 ! LCS_WINDOWS_COLOR_SPACE
! Remaining fields ignored: we will set them to zero.
   bmp_header%bV5EndPoints = CIEXYZTRIPLE( &
      CIEXYZ(0,0,0), &
      CIEXYZ(0,0,0), &
      CIEXYZ(0,0,0))
   bmp_header%bV5GammaRed    = 0
   bmp_header%bV5GammaGreen  = 0
   bmp_header%bV5GammaBlue   = 0
   bmp_header%bV5Intent      = 0
   bmp_header%bV5ProfileData = 0
   bmp_header%bV5ProfileSize = 0
   bmp_header%bV5Reserved    = 0
! *************************************
! Create palette
! Black & white 2-color palette
   allocate(palette(0:2**bmp_header%bV5BitCount-1))
! 0 = black
   palette(0) = RGBQUAD(0_ik1,0_ik1,0_ik1,0_ik1)
! 1 = white
   palette(1) = RGBQUAD(-1_ik1,-1_ik1,-1_ik1,0_ik1)

! *************************************
! Define image dimensions
! Calculate image width
   bmp_header%bV5Width = nint(width*DPI)
   RowLength = 4*((bmp_header%bV5Width*bmp_header%bV5BitCount+31)/32)
   bmp_header%bV5Width = RowLength*8/bmp_header%bV5BitCount
   RowLength = 4*((bmp_header%bV5Width*bmp_header%bV5BitCount+31)/32)
! Calculate image height
   bmp_header%bV5Height = nint(height*DPI)
   ColHeight = bmp_header%bV5Height
! Allocate image
   allocate(image(0:Rowlength-1,0:ColHeight-1))
! Store size of image in bytes
   bmp_header%bV5SizeImage = size(image)
! Store offset in file of DIB pixel bits
   bfOffsetBits = size(transfer(bmp_file_header,(/1_ik1/)))+ &
      bmp_header%bV5Size+ &
      size(transfer(palette,(/1_ik1/)))
!   bmp_file_header%bfOffsetBits = transfer(bfOffsetBits, &
!      bmp_file_header%bfOffsetBits)
! Work around ftn95 bug
   bmp_file_header%bfOffsetBits = transfer(bfOffsetBits, &
      (/1_ik2/))
   bfSize = bfOffsetBits+size(image)
!   bmp_file_header%bfSize = transfer(bfSize, &
!      bmp_file_header%bfSize)
! Work around ftn95 bug
   bmp_file_header%bfSize = transfer(bfSize, &
      (/1_ik2/))
! Work on image
!   lambda = 0.6563e-6_dp ! Red light
!lambda = 1.0e-5_dp
   write(*,'(a)',advance='no') 'Enter the wavelength in meters:> '
   read(*,*) lambda
   if(lambda <= 0) then
      write(flt_err,*) lambda
      write(*,'(3a)') 'Invalid wavelength: ',trim(adjustl(flt_err)),' meters'
      write(*,'(a)') 'Program aborting'
      stop
   end if
!   f = 0.3_dp ! Focal lenth = 30 cm
!f = 0.5_dp
   write(*,'(a)',advance='no') 'Enter the focal length in meters:> '
   read(*,*) f
   if(f <= 0) then
      write(flt_err,*) f
      write(*,'(3a)') 'Invalid focal length: ',trim(adjustl(flt_err)),' meters'
      write(*,'(a)') 'Program aborting'
      stop
   end if
   write(*,'(a)',advance='no') 'Use dithering? (default=no):> '
   read(*,'(a)',advance='no',eor=10,size=line_size) line
   read(*,'()')
10 continue
   if(line_size == 0) line(1:1) = 'N'
   line = adjustl(line)
   dither = scan(line(1:1),'tTyY') /= 0
   if(dither) then
      apodize = .FALSE.
   else
      write(*,'(a)',advance='no') 'Apodize? (default=no):> '
      read(*,'(a)',advance='no',eor=20,size=line_size) line
      read(*,'()')
20    continue
      if(line_size == 0) line(1:1) = 'N'
      line = adjustl(line)
      apodize = scan(line(1:1),'tTyY') /= 0
   end if
   r0 = f/sqrt((1/(PelsPerMeter*lambda))**2+1)
   sigma = r0
! The following tests for g95 and others
   if(is_ftn95()) then ! Use ftn95-style stream I/O
      access = 'transparent'
      form = 'unformatted'
   else if(is_g95()) then ! Use f03-compatible stream I/O
      access = 'stream'
      form = 'unformatted'
   else ! Use values appropriate for lf95 express 7.10.02,
        ! ifort Package ID: W_FC_C_9.0.029, or cvf 6.6c3
      access = 'sequential'
      form = 'binary'
   end if
   open(10,file='fire.bmp',access=access,form=form,status='replace')
   write(10) bmp_file_header
   write(10) bmp_header
   write(10) palette
   x0 = bmp_header%bV5Width/2
   y0 = bmp_header%bV5Height/2
! Initialize image to white
   image = -1_ik1
   if(dither) then
      call random_seed()
      do i = 0, bmp_header%bV5Width-1
         do j = 0, bmp_header%bV5Height-1
            call random_number(harvest)
            r = sqrt(((i-x0+harvest(1))/real(PelsPerMeter,dp))**2+ &
                     ((j-y0+harvest(2))/real(PelsPerMeter,dp))**2+ &
                     f**2)
            if(abs(mod(r/lambda,1.0_dp)-0.5_dp) < 0.25_dp) then
               image(i/8,j) = iand(image(i/8,j),not(ishft(1_ik1,7-mod(i,8))))
            else
               image(i/8,j) = ior(image(i/8,j),ishft(1_ik1,7-mod(i,8)))
            end if
         end do
      end do
   else if(apodize) then
      do i = 0, bmp_header%bV5Width-1
         do j = 0, bmp_header%bV5Height-1
            r12 = ((i-x0)/real(PelsPerMeter,dp))**2+ &
                 ((j-y0)/real(PelsPerMeter,dp))**2
            r = sqrt(r12+f**2)
            r1 = sqrt(r12)
            if(abs(mod(r/lambda,1.0_dp)-0.5_dp) < 0.25_dp*exp(-(r1/sigma)**2)) then
               image(i/8,j) = iand(image(i/8,j),not(ishft(1_ik1,7-mod(i,8))))
            else
               image(i/8,j) = ior(image(i/8,j),ishft(1_ik1,7-mod(i,8)))
            end if
         end do
      end do
   else
      do i = 0, bmp_header%bV5Width-1
         do j = 0, bmp_header%bV5Height-1
            r12 = ((i-x0)/real(PelsPerMeter,dp))**2+ &
                 ((j-y0)/real(PelsPerMeter,dp))**2
            r = sqrt(r12+f**2)
            r1 = sqrt(r12)
!            if(abs(mod(r/lambda,1.0_dp)-0.5_dp) < 0.25_dp .OR. &
!               abs(r1-r0) < 3.0_dp/PelsPerMeter) then ! Plot circle
            if(abs(mod(r/lambda,1.0_dp)-0.5_dp) < 0.25_dp) then
               image(i/8,j) = iand(image(i/8,j),not(ishft(1_ik1,7-mod(i,8))))
            else
               image(i/8,j) = ior(image(i/8,j),ishft(1_ik1,7-mod(i,8)))
            end if
         end do
      end do
   end if
   write(10) image
   close(10)
end program fire
! End of file: fire2.f90