ラプラス方程式　アニメーション gif 化

#能書き
以前、二次元ラプラス方程式をキャラクタ・グラフィックスで表示してみましたが、それを内容そのままで出力形式をアニメーション gif にしてみます。

https://qiita.com/cure_honey/items/b8b106a309e505b654f6

アニメーション gif のルーチンは、 fortran66 さんから借りてくることにします。初期化と終了処理を除けば、以前のキャラクター出力の所を少しいじるだけで直せます。色の設定が難しいですね。

http://fortran66.hatenablog.com/entry/2018/05/09/014127

ラプラス方程式の解は、微分方程式を反復法で求めています。ここで mask 付きの forall 構文を用いることで、微分方程式の反復解法の記述と境界条件の記述とを分離することができます。

プログラムでは往々にして、解法の記述と境界条件の記述とが渾然一体となってしまうので、この方式だとスッキリ気分を味わえます。
##出力結果
外側のパイプが 100 度、内側のパイプが 0 度、パイプ間の媒質が 20 度という初期値で計算しています。図では温度を 10 度刻みで色の濃淡で表わしています。
###同心円の場合
```fortran
  ! define shape : concentric circle  
      mask = .false.       
      forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mask(i, j) = .true.
      forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 <  50**2) mask(i, j) = .false.
  ! initial value
      mesh = 100.0
      forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mesh(i, j) = 20.0    
      forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 <  50**2) mesh(i, j) =  0.0    
```
![heat-con.gif](https://qiita-image-store.s3.amazonaws.com/0/44522/c618acac-7fa9-827e-e950-2cda8a67cf61.gif)

###中心のずれた円の場合
```fortran
      ! define shape : eccentric pipe   
          mask = .false.       
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mask(i, j) = .true.
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 <  50**2) mask(i, j) = .false.
      
      ! initial value
          mesh = 100.0
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mesh(i, j) = 20.0    
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 <  50**2) mesh(i, j) =  0.0
```
![heat-ecc.gif](https://qiita-image-store.s3.amazonaws.com/0/44522/756f9ba6-b184-f08a-7430-d4524a94a69c.gif)
###内側のパイプが矩形の場合
```fortran
  ! define shape : circle and square   
      mask = .false.       
      forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2     ) mask(i, j) = .true.
      forall (i = 1:n, j = 1:n, abs(i - 351) < 75 .and. abs(j - 351) < 75) mask(i, j) = .false.
  ! initial value
      mesh = 100.0
      forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2     ) mesh(i, j) = 20.0    
      forall (i = 1:n, j = 1:n, abs(i - 351) < 75 .and. abs(j - 351) < 75) mesh(i, j) =  0.0
```
![heat-sq.gif](https://qiita-image-store.s3.amazonaws.com/0/44522/4d0de9a5-33ad-6199-7765-3de60964d361.gif)

##プログラム・リスト
プログラムでは、Fortran2008 の機能を使っています。forall(integer:: ) の構文で局所変数を使っているので、Intel Fortran の Ver.19 beta 版でないと動きません。（integer を普通に変数宣言すれば F95 形式になります。)なお、コンパイラのオプションで、semantics を Fortran2003 対応にしておかないと、微妙なメモリー・エラーが出たり出なかったりします。

Fortran2008 は block 構文が導入されて、局所変数が使えるだけでなく、use 文も使えるようになりました。サブルーチン化するほどでもないが局所的にまとめておきたい、たとえば今回の様なグラフィックスのような目的に便利な機能です。

```fortran
    !
    ! animation gif modules borrowed from http://fortran66.hatenablog.com/entry/2018/05/09/014127
    !
    module m_gif_types
        use, intrinsic :: iso_fortran_env
        implicit none    
        private
        public :: t_gif_header, t_image_block, t_graphic_control_extension, t_application_extension 
      
      ! constructors (interface)
        interface t_gif_header
            procedure :: init_gif_header
        end interface 

        
        interface t_image_block
            procedure :: init_image_block
        end interface 
      
        
        interface t_graphic_control_extension
            procedure :: init_graphic_control_extension 
        end interface
      
        
        interface t_application_extension
            procedure :: init_application_extension
        end interface  
      !
        type :: t_gif_header
            sequence  
            character(3)   :: signature               = 'GIF'
            character(3)   :: version                 = '89a' ! '87a'
            integer(int16) :: width     
            integer(int16) :: height    
            integer(int8)  :: pck                     = int(B'10001010', int8) !1:Global Col. 3:Color Res. 1:Sort Flg. 3:Size of Global Col.
            integer(int8)  :: background_color_index  = 0
            integer(int8)  :: pixel_aspect_ratio      = 0
        end type  

        
        type :: t_image_block
            sequence
            integer(int8)  :: image_separator         = int(Z'2C', int8)
            integer(int16) :: image_left_position     = 00
            integer(int16) :: image_top_position      = 00
            integer(int16) :: image_width  
            integer(int16) :: image_height 
            integer(int8)  :: pck                     = int(B'00000000', int8)
            ! local color table
            !integer(int8) :: LZW_minimum_code_size 
            !integer(int8) :: block_size  
            !integer(int8) :: image_data (:)
            !integer(int8) :: block_terminator        = int(Z'00', int8)  
        end type t_image_block

        
        type :: t_graphic_control_extension
            sequence  
            integer(int8)  :: extention_introducer    = int(Z'21', int8)
            integer(int8)  :: graphic_control_label   = int(Z'F9', int8)
            integer(int8)  :: lock_size               = int(Z'04', int8)
            integer(int8)  :: pck                     = int(Z'00', int8)
            integer(int16) :: delay_time              = 50 ! msec
            integer(int8)  :: transparent_color_index = 0
            integer(int8)  :: block_terminator        = int(Z'00', int8)
        end type t_graphic_control_extension
      
        
        type :: t_application_extension ! for animation gif
            sequence  
            integer(int8)  :: extension_intrducer     = int(Z'21', int8)
            integer(int8)  :: extension_label         = int(Z'FF', int8)
            integer(int8)  :: block_size_01           = int(Z'0B', int8)
            character(len = 8) :: application_identifier          = 'NETSCAPE'
            character(len = 3) :: application_authentication_code = '2.0'
            integer(int8)  :: block_size_02           = int(Z'03', int8) ! 0:block terminator
            integer(int8)  :: n                       = int(Z'01', int8) !application_data
            integer(int16) :: nloop                   = 0                ! 0:unlimited loop 
            integer(int8)  :: block_terminator        = int(Z'00', int8)
        end type t_application_extension   
    
    contains
    
        ! constructors (implementation)
        pure type(t_gif_header) function init_gif_header(nx, ny, ngcol) result(res)
            integer, intent(in) :: nx, ny, ngcol
            res%width  = nx
            res%height = ny
            res%pck = iand(res%pck, int(B'11111000', int8)) + int(mod(ngcol, 8), int8) !least 3bits:size of global color:2**(ngcol + 1) 2..256
        end function init_gif_header
        
        
        pure type(t_image_block) function init_image_block(nx, ny) result(res)
            integer, intent(in) :: nx, ny
            res%image_width  = nx
            res%image_height = ny
        end function init_image_block
        
         
        pure type(t_graphic_control_extension) function init_graphic_control_extension(it) result(res)
            integer, intent(in) :: it
            res%delay_time = int(it, int16) ! msec 
        end function  init_graphic_control_extension
        
    
        pure type(t_application_extension) function init_application_extension(nloop) result(res)
            integer, intent(in) :: nloop
            res%nloop = int(nloop, int16)
        end function init_application_extension  
    end module m_gif_types 

        
    module m_gif_lzw
        use, intrinsic :: iso_fortran_env
        use :: m_gif_types
        implicit none
        private
        public :: encoder
        type :: t_enc
            integer :: kbits = 0, id = 0
        end type
    contains
    
        subroutine enc(irgb, nbits, code)                   ! gif LZW 
            integer, intent(in) :: irgb(:), nbits
            type (t_enc), intent(out) :: code(:)
            integer :: dict(0:2**12 - 1) ! 0:4095
            integer :: kbits, kd, m0                        ! <== m0 state variable
            integer :: ip, ic, ienc
            call clear_dict()                               ! clear dictionary 
            ic = 1
            code(ic) = t_enc(kbits, 2**nbits)               ! clear_code 
            m0 = irgb(1)
            do ip = 2, size(irgb)
                call subenc(ienc, irgb(ip)) 
                if (ienc /= -1) then                        ! new word 
                    ic = ic + 1
                    code(ic) = t_enc(kbits, ienc)
                end if  
                if (kd == 2**12 - 1) then                   ! dictionary full
                    ic = ic + 1  
                    code(ic) = t_enc(kbits, 2**nbits)       ! clear_code
                    call clear_dict()                       ! clear dictionary 
                end if
            end do
            ic = ic + 1
            code(ic) = t_enc(kbits, m0)
            ic = ic + 1 
            code(ic) = t_enc(kbits, dict(2**nbits + 1))     ! end_code
            return
        contains
        
            subroutine subenc(ienc, m1)
                integer, intent(out) :: ienc
                integer, intent(in ) :: m1
                integer :: k, id
                k = ishft(m0 + 1, 16) + m1                  ! m0 + 1 to avoid 00...0 degeneracy 
                id = findloc(dict, k, dim = 1)              ! dictionary is 0-based 
                if (id == 0) then                           ! not found in the dictionary
                    kd = kd + 1
                    dict(kd) = k
                    if (kd == 2**kbits + 1) kbits = kbits + 1
                    ienc = m0
                    m0 = m1
                else                                        ! found in the dictionary
                    ienc = -1
                    m0 = id - 1                             ! because dictionary is 0-based, m0 must be shifted by 1
                end if    
              end subroutine subenc
             
              subroutine clear_dict()
                  kbits = nbits + 1
                  kd = 2**nbits + 1 ! dictinary defined
                  dict = 0
                  forall(integer :: i = 0:kd) dict(i) = i   ! dict(nbit):clear code; dict(nbit+1):end code 
              end subroutine clear_dict
          end subroutine enc
      
        
          subroutine encoder(irgb, nbits, icode)
              integer, intent(in ) :: irgb(:), nbits
              integer(int8), allocatable, intent(out) :: icode(:)
              type(t_enc), allocatable :: code(:)
              integer :: i, j, k, ib, ic, nb  
              allocate(code(size(irgb) + 100))              ! ? large enough ?
              call enc(irgb, nbits, code)                   ! compress color code to LZW code
              nb = ceiling(sum(code(:)%kbits) / 8.0)        !   required bytes  
              allocate(icode(nb))
              icode = 0                                     ! pack LZW code to bit string    
              k = 0
              do i = 1, size(code)                          
                  do j = 1, code(i)%kbits   
                      if (btest(code(i)%id, j - 1)) then 
                      ic = k / 8 + 1 
                      ib = mod(k, 8)
                      icode(ic) = ibset(icode(ic), ib)
                  end if
                      k = k + 1
                  end do  
              end do            
          end subroutine encoder 
      end module m_gif_lzw
     
    
      module m_anime_gif
          use, intrinsic :: iso_fortran_env
          use :: m_gif_types
          use :: m_gif_lzw
          implicit none
          private 
          public :: t_gif
          
          interface t_gif 
              procedure :: init_gif
          end interface
          
          type :: t_gif
              integer :: iw, nbits
              integer, allocatable :: irgb2(:, :)
          contains
              procedure :: wr_gif_img, gif_close
          end type t_gif
          
      contains
    
          impure type(t_gif) function init_gif(fn, nx, ny, n_global_color, global_color_table, loop) result(this)
              character(*), intent(in) :: fn
              integer, intent(in) :: nx, ny, n_global_color, global_color_table(:)
              integer, intent(in), optional :: loop
              type(t_gif_header) :: head
              type(t_application_extension) :: anime
              integer :: nloop = 0
              if (present(loop)) nloop = loop
              associate(iw => this%iw, nbits => this%nbits)
                  nbits = max(2, n_global_color + 1) ! 1..8 -> 2,2,3,4,5,6,7,8
                  allocate(this%irgb2(nx, ny))
                  head  = t_gif_header(nx, ny, n_global_color)
                  anime = t_application_extension(nloop)
                  open(newunit = iw, file = fn, access = 'stream') 
                  write(iw) head      
                  write(iw) int(global_color_table, int8)
                  write(iw) anime
              end associate  
          end function init_gif

          subroutine wr_gif_img(this, itime)
              class(t_gif), target, intent(in) :: this
              integer, intent(in) :: itime
              type(t_graphic_control_extension), allocatable :: ext
              type(t_image_block), allocatable :: img
              integer(int8), allocatable :: icode(:)
              integer :: i, nbyte  
              integer, pointer :: irgb(:)
              associate(iw => this%iw, nbits => this%nbits, irgb2 => this%irgb2)
                  ext = t_graphic_control_extension(itime)
                  img = t_image_block(size(irgb2, 1), size(irgb2, 2))
                  irgb(1:size(irgb2)) => irgb2      ! 1D => 2D 
                  call encoder(irgb, nbits, icode)  ! color code to GIF LZW code
                  write(iw) ext                     ! write bit strings to file
                  write(iw) img     
                  write(iw) int(nbits, int8)        ! LZW_minimum_code_size = 2,2,3,4,5,6,7,8bits 
                  do i = 1, size(icode), 255         
                    nbyte = min(i + 255, size(icode)) - i
                    write(iw) achar(nbyte)             
                    write(iw) icode(i:i+nbyte-1)
                  end do  
                  write(iw) achar(00)               ! block_terminator
              end associate
          end subroutine wr_gif_img
      
    
          subroutine gif_close(this)
              class(t_gif), intent(in out) :: this
              associate(iw => this%iw)
                  write(iw) achar(int(Z'3B'))                       
                  close(iw)
              end associate  
          end subroutine gif_close

      end module m_anime_gif
    
    
      program Heat
          implicit none
          integer, parameter :: n = 501
          real    :: mesh(0:n + 1, 0:n + 1)
          logical :: mask(n, n)
          integer :: iter
          
      ! define shape : eccentric pipe   
          mask = .false.       
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mask(i, j) = .true.
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 <  50**2) mask(i, j) = .false.
      
      ! initial value
          mesh = 100.0
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mesh(i, j) = 20.0    
          forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 <  50**2) mesh(i, j) =  0.0
      
      ! solve Laplace's equation by iterative method
          block 
              use m_anime_gif
              type(t_gif), allocatable :: fig
              integer :: n_global_color               ! 0..7 ! 2..256 colors (2^(n_global_color + 1)) 
              integer, allocatable :: global_color_table(:)
           
              n_global_color = 7
              allocate(global_color_table(3 * 2**(n_global_color + 1)))
              global_color_table = 0
              forall (integer :: i = 0:254) global_color_table(3*i+1:3*i+3) = mod(2 * i, 255)
              fig = t_gif('heat.gif', n+2, n+2, n_global_color, global_color_table)
              !
              do iter = 0, 100000
                  if (mod(iter, 5000) == 0) then        ! gif animation
                      fig%irgb2 = 10 * int(mesh / 10.0)
                      call fig%wr_gif_img(20)
                  end if  
                  forall (integer:: i = 1:n, j = 1:n, mask(i, j)) & 
                      mesh(i, j) = (mesh(i - 1, j) + mesh(i + 1, j) + mesh(i, j - 1) + mesh(i, j + 1)) / 4
              end do
              call fig%gif_close()
              deallocate(fig)
          end block  
      end program Heat
```