module para
  implicit none
  integer :: i
  integer,parameter :: N = 1001
  real(8) :: dx,pi
end module para
!*********************************************************************
program euler
  use para
  implicit none
  real(8),dimension(N) :: x,y,f
  
  call calc_param(x,y)
  call calc_Lhand(x,f)
  call calc_y(x,y,f)
  call init
  call viewport(0.2,0.2,0.8,0.8)
  call xyworld(0.0,-1.2,2.4,1.2)
  call frame
  call plotd(x,y)
  call fin
  stop
  
end program euler
!*********************************************************************
subroutine calc_param(x,y)
  use para
  implicit none
  real(8),dimension(N) :: x,y
  
  pi = 4.d0*datan(1.d0)
  dx = 2.d0*pi/dble(N-1)
  
  do i = 1,N
    x(i) = dble(i-1)*dx
  end do
  
!////inicialize & boundary condition
  y(:) = 0.d0
  y(1) = 1.d0
  
  return
end subroutine calc_param
!*********************************************************************
subroutine calc_Lhand(x,f)
  use para
  implicit none
  real(8),dimension(N) :: x,f
  
  do i = 1,N
    f(i) = -dsin(x(i))
  end do
  
  return
end subroutine calc_Lhand
!*********************************************************************
subroutine calc_y(x,y,f)
  use para
  implicit none
  real(8),dimension(N) :: x,y,f
  
  do i = 2,N
    y(i) = y(i-1)+dx*f(i-1)
  end do
  
  do i = 1,N
    write(*,*) 'x=',x(i),'y=',y(i)
  end do
  
  return
end subroutine calc_y
!*********************************************************************
subroutine frame
  use para
  implicit none
  
  call line1(0.0,0.0,2.0,0.0)
  call line1(real(pi/4.d0),-0.02,real(pi/4.d0),0.02)
  call line1(real(pi/2.d0),-0.02,real(pi/2.d0),0.02)
  call line1(0.0,0.5,0.04,0.5)
  call line1(0.0,-0.5,0.04,-0.5)
  call texty(0.0,0.0,1,'y')
  call texty(0.0,1.0,3,'1.0')
  call texty(0.0,-1.0,4,'-1.0')
  call texty(0.0,0.5,3,'0.5')
  call texty(0.0,-0.5,4,'-0.5')
  call text(1.0,-1.1,1,'x')
  call rect1(0.0,-1.0,2.0,1.0)

  call setchar(1,15)
  call textx(real(pi/4.d0),-0.03,4,'pi/4')
  call textx(real(pi/2.d0),-0.03,4,'pi/2')
  call stroke
  
  return
end subroutine frame
!*********************************************************************
subroutine plotd(x,y)
  use para
  implicit none
  integer :: ip = 1
  real(8),dimension(N) :: x,y
  real(4),dimension(N) :: x0,y0
  
  x0 = real(x)
  y0 = real(y)
  
  do i = 1,N
    if(x0(i) > 2.e0) cycle
    call plot(x0(i),y0(i),ip+2)
    ip = 0
  end do
  
  call stroke
  
  return
  
end subroutine plotd