計算は手抜きのオイラー法で。話を散乱対を含む二次元面内に限定します。
実行結果
中央付近の四角い枠の中心に固定された正の点電荷があるとき、左側から縦の位置を乱数で決めた、右に水平に同じ速度を持つ荷電粒子を同時に発射したときの散乱のようすです。射ち込まれる荷電粒子間の相互作用は考えません。射ち込まれる荷電粒子の電荷が正の場合と負の場合の計算結果を示しています。
１００個の粒子を同時に、１ステップづつ進ませています。ELEMENTAL なサブルーチンを使うことで、明示的に個数に関するループを回さなくて済みます。
粒子の個数は宣言文で、粒子位置の更新は配列名のみで個数を非明示に行います。
type(t_particle) :: p(100), q(100) ... call p%next(0.1) call q%next(0.1)
絵を描いている途中でもコンソールからの入出力は可能です。
ソース・プログラム
module m_oop implicit none type :: t_rgb integer :: ir, ig, ib end type t_rgb type, abstract :: t_device character(len = 80) :: title = 'Plotter' integer :: nsize_x = 640, nsize_y = 480 integer :: line_width = 1 type (t_rgb) :: rgb = t_rgb(0, 0, 0) contains procedure (device_on), deferred, pass :: on procedure (device_off), deferred, pass :: off procedure (device_show), deferred, pass :: show procedure (device_pen), deferred, pass :: pen procedure (device_lineTo), deferred, pass :: lineTo procedure (device_moveTo), deferred, pass :: moveTo procedure (device_dot), deferred, pass :: dot end type t_device abstract interface subroutine device_on(self) import :: t_device class(t_device), intent(in out) :: self end subroutine device_on subroutine device_off(self) import :: t_device class(t_device), intent(in) :: self end subroutine device_off subroutine device_show(self) import :: t_device class(t_device), intent(in) :: self end subroutine device_show subroutine device_pen(self, line_width, rgb) import :: t_device, t_rgb class(t_device), intent(in out) :: self integer, intent(in), optional :: line_width type (t_rgb), intent(in), optional :: rgb end subroutine device_pen subroutine device_lineTo(self, ix, iy) import :: t_device class(t_device), intent(in) :: self integer, intent(in) :: ix, iy end subroutine device_lineTo subroutine device_moveTo(self, ix, iy) import :: t_device class(t_device), intent(in) :: self integer, intent(in) :: ix, iy end subroutine device_moveTo subroutine device_dot(self, ix, iy, icol) import :: t_device, t_rgb class(t_device), intent(in) :: self integer, intent(in) :: ix, iy integer, intent(in) :: icol ! type (t_rgb), intent(in) :: rgb end subroutine device_dot end interface end module m_oop module m_win32 use ifwina use ifwinty use ifmt, only : RTL_CRITICAL_SECTION use m_oop implicit none type :: t_wnd integer (HANDLE) :: hWnd = 0 integer (HANDLE) :: hDC = 0 integer (LPINT) :: hThread = 0 integer (LPDWORD):: id = 0 integer (HANDLE) :: hPen = 0 type (RTL_CRITICAL_SECTION) :: lpCriticalSection = RTL_CRITICAL_SECTION(0,0,0,0,0,0) end type t_wnd type, extends(t_device) :: t_win32 type (t_wnd) :: wnd contains procedure, pass :: on => gr_on procedure, pass :: off => gr_off procedure, pass :: show => gr_show procedure, pass :: pen => gr_pen procedure, pass :: lineTo => gr_lineTo procedure, pass :: moveTo => gr_moveTo procedure, pass :: dot => gr_dot end type t_win32 integer (DWORD), save :: iTls = 0 integer, save :: nwin = 0 contains !-------------------------------------------------------------------------------- integer(4) function WinMain( hInstance, nCmdShow, win32 ) implicit none integer (HANDLE), intent(in) :: hInstance integer (SINT) , intent(in) :: nCmdShow type (t_win32) , intent(in) :: win32 type (T_WNDCLASS) :: wc type (T_MSG) :: mesg integer (HANDLE) :: hWndMain integer (BOOL) :: iretb character (LEN = 256) :: ClassName = 'Fortran'//char(0) integer :: iwindow_frame_x, iwindow_frame_y logical, save :: first = .true. ! Init Main window iwindow_frame_x = 2 * GetSystemMetrics(SM_CXFIXEDFRAME) !side line = 6, title bar = 25 iwindow_frame_y = 2 * GetSystemMetrics(SM_CYFIXEDFRAME) + GetSystemMetrics(SM_CYCAPTION) ! if (first) then WinMain = -1 ! Error code wc%lpszClassName = loc(ClassName) ! non-standard Fortran :: LOC(xxx) = TRANSFER(C_LOC(xxx), iii) wc%lpfnWndProc = loc(MainWndProc) ! CALLBACK procedure name wc%style = ior(CS_VREDRAW , CS_HREDRAW) wc%hInstance = hInstance wc%hIcon = NULL wc%hCursor = LoadCursor( NULL, IDC_ARROW ) wc%hbrBackground = ( COLOR_WINDOW + 1 ) if ( RegisterClass(wc) == 0 ) return ! initialize window first = .false. end if ! Init instance WinMain = -2 ! Error code hWndMain = CreateWindowEx( 0, ClassName, & trim(win32%title)//char(0), & int(ior(WS_OVERLAPPED, WS_SYSMENU)), & CW_USEDEFAULT, CW_USEDEFAULT, & win32%nsize_x + iwindow_frame_x, & win32%nsize_y + iwindow_frame_y, & 0, 0, & hInstance, & LOC(win32%wnd) ) ! non-standard if (hWndMain == 0) return iretb = ShowWindow( hWndMain, nCmdShow ) iretb = UpdateWindow( hWndMain ) ! Message Loop do while ( GetMessage (mesg, NULL, 0, 0) ) iretb = TranslateMessage( mesg ) iretb = DispatchMessage( mesg ) end do WinMain = mesg%wParam return end function WinMain !------------------------------------------------------------------------------------- integer (LRESULT) function MainWndProc( hWnd, mesg, wParam, lParam ) !DEC$ ATTRIBUTES STDcall, DECORATE, ALIAS : 'MainWndProc' :: MainWndProc integer (HANDLE) , intent(in) :: hWnd integer (UINT) , intent(in) :: mesg integer (fwParam), intent(in) :: wParam integer (flParam), intent(in) :: lParam ! integer (HANDLE) :: hDC, hBmp integer (BOOL) :: iretb type (T_PAINTSTRUCT) :: ps type (T_RECT) :: rc ! TYPE (T_CREATESTRUCT) :: cs TYPE (t_wnd) :: wnd POINTER (p_cs, cs) ! non-standard POINTER (p_wnd, wnd) ! non-standard ! MainWndProc = 0 select case ( mesg ) case (WM_CREATE) p_cs = lParam p_wnd = cs%lpCreateParams iretb = TlsSetValue(iTls, p_wnd) wnd%hWnd = hWnd hDC = GetDC(hWnd) wnd%hDC = CreateCompatibleDC(hDC) iretb = GetClientRect(hWnd, rc) hBmp = CreateCompatibleBitmap(hDC, rc%right - rc%left, rc%bottom - rc%top) iretb = SelectObject(wnd%hDC, hBmp) iretb = PatBlt(wnd%hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, WHITENESS) iretb = ReleaseDC(hWnd, hDC) iretb = DeleteObject(hBmp) case (WM_DESTROY) p_wnd = TlsGetValue(iTls) call EnterCriticalSection( loc(wnd%lpCriticalSection) ) iretb = DeleteObject( wnd%hDC ) call PostQuitMessage( 0 ) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) case (WM_PAINT) p_wnd = TlsGetValue(iTls) call EnterCriticalSection( loc(wnd%lpCriticalSection) ) hDC = BeginPaint( wnd%hWnd, ps ) iretb = GetClientRect( wnd%hWnd, rc ) iretb = BitBlt(hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, wnd%hDC, 0, 0, SRCCOPY) iretb = endPaint( wnd%hWnd, ps ) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) case (WM_RBUTTONUP) p_wnd = TlsGetValue(iTls) call EnterCriticalSection( loc(wnd%lpCriticalSection) ) iretb = DeleteObject( wnd%hDC ) call PostQuitMessage( 0 ) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) case default MainWndProc = DefWindowProc( hWnd, mesg, wParam, lParam ) end select return end function MainWndProc !------------------------------------------------------------------------------------- subroutine gr_on(self) use IFMT, only : CreateThread ! multithread module class(t_win32), intent(in out) :: self integer (BOOL) :: iretb integer (HANDLE) :: hBmp type (T_RECT) :: rc associate(wnd => self%wnd) call InitializeCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC if (nwin == 0) iTls = TlsAlloc() nwin = nwin + 1 wnd%hThread = CreateThread(NULL, 0, Thread_Proc, NULL, CREATE_SUSPENDED, wnd%id) iretb = SetThreadPriority(wnd%hThread, THREAD_PRIORITY_BELOW_NORMAL) iretb = ResumeThread(wnd%hThread) call sleep(100) ! wait for Window initialization iretb = GetClientRect(wnd%hWnd, rc) hBmp = CreateCompatibleBitmap(wnd%hDC, rc%right - rc%left, rc%bottom - rc%top) iretb = SelectObject(wnd%hDC, hBmp) iretb = DeleteObject(hBmp) iretb = PatBlt(wnd%hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, WHITENESS) wnd%hPen = CreatePen(PS_SOLID, 1, 0) end associate return contains integer (LONG) function Thread_Proc(lp_ThreadParameter) !DEC$ ATTRIBUTES STDcall, ALIAS:"_thread_proc" :: Thread_Proc integer (LPINT), intent(in) :: lp_ThreadParameter integer :: hInst hInst = GetModuleHandle(NULL) Thread_Proc = WinMain(hInst, SW_SHOWNORMAL, self) return end function Thread_Proc end subroutine gr_on !------------------------------------------------------------------------------------- subroutine gr_off(self) class(t_win32), intent(in) :: self integer (BOOL) :: iretb integer (DWORD) :: iwait associate(wnd => self%wnd) iwait = INFINITE iretb = InvalidateRect(wnd%hWnd, NULL, FALSE) iretb = DeleteObject(wnd%hPen) iretb = WaitForSingleObject(wnd%hThread, iwait) iretb = CloseHandle(wnd%hThread) iretb = PostMessage(wnd%hWnd, WM_DESTROY, NULL, NULL) ! wnd%hThread = NULL call DeleteCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC end associate nwin = nwin - 1 if (nwin ==0) iretb = TlsFree(iTls) return end subroutine gr_off !------------------------------------------------------------------------------------- subroutine gr_show(self) class(t_win32), intent(in) :: self integer (BOOL):: iretb associate(wnd => self%wnd) call EnterCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC iretb = InvalidateRect(wnd%hWnd, NULL, FALSE) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC end associate return end subroutine gr_show !------------------------------------------------------------------------------------- subroutine gr_pen(self, line_width, rgb) class(t_win32), intent(in out) :: self integer, intent(in), optional :: line_width type (t_rgb), intent(in), optional :: rgb integer (BOOL) :: iretb associate( rgb_ => self%rgb, line_width_ => self%line_width, wnd => self%wnd ) if ( present(rgb) ) rgb_ = rgb if ( present(line_width) ) line_width_ = line_width call EnterCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC iretb = DeleteObject(wnd%hPen) wnd%hPen = CreatePen(PS_SOLID, line_width_, irgb(rgb_)) iretb = SelectObject(wnd%hDC, wnd%hPen) iretb = MoveToEx(wnd%hDC, 0, 0, NULL) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC end associate return contains integer function irgb(rgb) type(t_rgb), intent(in) :: rgb irgb = rgb%ir + (rgb%ig + (rgb%ib * 256)) * 256 return end function irgb end subroutine gr_pen !---------------------------------------------------------------- subroutine gr_moveTo(self, ix, iy) class(t_win32), intent(in) :: self integer, intent(in) :: ix, iy integer (BOOL):: iretb associate(wnd => self%wnd) call EnterCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC iretb = MoveToEx(wnd%hDC, ix, iy, NULL) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC end associate return end subroutine gr_moveTo !---------------------------------------------------------------- subroutine gr_lineTo(self, ix, iy) class(t_win32), intent(in) :: self integer, intent(in) :: ix, iy integer (BOOL):: iretb associate(wnd => self%wnd) call EnterCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC iretb = LineTo(wnd%hDC, ix, iy) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC end associate return end subroutine gr_lineTo !------------------------------------------------------------------------------------- subroutine gr_dot(self, ix, iy, icol) class(t_win32), intent(in) :: self integer, intent(in) :: ix, iy, icol integer (BOOL):: iretb associate(wnd => self%wnd) call EnterCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC iretb = SetPixel(wnd%hDC, ix, iy, icol) call LeaveCriticalSection( loc(wnd%lpCriticalSection) ) ! non-standard Fortran :: LOC end associate return end subroutine gr_dot !---------------------------------------------------------------- end module m_win32 module m_plot use m_oop use m_win32 implicit none private public :: t_rgb, t_device, t_win32, t_wnd end module m_plot module m_particle implicit none type :: t_particle real :: r(3) = 0.0 real :: v(3) = 0.0 real :: xmass = 1.0, q = +1.0 contains procedure :: next procedure :: force end type t_particle contains pure elemental subroutine next(this, dt) class (t_particle), intent(in out) :: this real, intent(in) :: dt this%v = this%v + dt * this%force() / this%xmass this%r = this%r + dt * this%v return end subroutine next pure function force(this) class (t_particle), intent(in) :: this real :: force(3) real :: r r = sqrt( sum(this%r**2) ) force = -this%q * this%r / r**3 return end function force end module m_particle program test use m_particle use m_plot implicit none type(t_particle) :: p(100), q(100) integer :: i, it, ix, iy class(t_device), allocatable ::fig1, fig2 call random_number(p%r(2)) call random_number(q%r(2)) p%r(1) = -2 p%v(1) = 1.0 q%r(1) = -2 q%v(1) = 1.0 q%q = -1.0 p%r(1) = 10.0 * (p%r(1) - 0.5) p%r(2) = 30.0 * (p%r(2) - 0.5) p%r(3) = 0.0 q%r(1) = 10.0 * (q%r(1) - 0.5) q%r(2) = 30.0 * (q%r(2) - 0.5) q%r(3) = 0.0 allocate(t_win32::fig1) fig1%title = 'Positive charge particles' fig1%nsize_x = 1000 fig1%nsize_y = 600 allocate(t_win32::fig2) fig2%title = 'Negative charge particles' fig2%nsize_x = 1000 fig2%nsize_y = 600 call fig1%on() call fig2%on() call square(fig1) call square(fig2) do it = 1, 700 print *, 'time =', it do i = 1, 100 ix = int(10 * p(i)%r(1) + 320) iy = int(10 * p(i)%r(2) + 300) call fig1%dot(ix, iy, irgb(0, 0, 255)) ix = int(10 * q(i)%r(1) + 320) iy = int(10 * q(i)%r(2) + 300) call fig2%dot(ix, iy, irgb(255, 0, 0)) end do call fig1%show() call fig2%show() call p%next(0.1) call q%next(0.1) call sleepqq(90) end do call fig1%off() call fig2%off() stop contains integer function irgb(ir, ig, ib) integer, intent(in) :: ir, ig, ib irgb = ir + 256 * (ig + 256 * ib) return end function irgb subroutine square(fig) class(t_device), intent(in) :: fig integer, parameter :: icx = 320, icy = 300, m = 3 call fig%moveto(icx + m, icy + m) call fig%lineto(icx + m, icy - m) call fig%lineto(icx - m, icy - m) call fig%lineto(icx - m, icy + m) call fig%lineto(icx + m, icy + m) return end subroutine square end program test
