program prot
!
!   Use BLAS subroutines DROTG and DROT to generate a rotation and calculate the 
!      rotation of a squared figure 
!
   implicit none
   integer, allocatable, dimension(:,:) :: plane
   real(8), allocatable, dimension(:,:,:) ::  square
   real(8), parameter :: PI=4.0D0*ATAN(1.0D0)
   real(8), dimension(100) :: x, y
   integer :: ld, n, lds, ns, rc, i, j, ix, iy, ii, jj
   real(8) :: da, db, c, s
   external :: drot
   
   n = 0
   print*,"Plane dimension?"
   read(*,*) n
   if ( n < 200 .OR. n > 1000 ) n = 1000
   ld = n
   allocate(plane(ld,n),stat=rc)
   if (rc /= 0 ) stop

   ns = n/4; lds = ns
   allocate(square(lds,ns,2),stat=rc)
   if (rc /= 0 ) stop
   
   plane = 0	! Initialize plane
   !
   ! Define squared figure ...
   !
   do i = 1, ns
      do j = 1, ns
      	if ( i == 1 .or. i == ns .or. j == 1 .or. j == ns ) then
      		square(i,j,1) = -ns/2 + i
      		square(i,j,2) = -ns/2 + j
      	end if
      end do
   end do
   	   
   ! ...  and trace it in plane
   call DrawFigure(ld,n,plane,lds,ns,square,0.5d0)
   
   print*,"Define rotation: da, db ? "
   read(*,*) da,db
!
   call drotg(da,db,c,s)	! get parameters C, S
   call drot(ns*ns,square(1,1,1),1,square(1,1,2),1,c,s)
!       
!	Draw new figure in plane
   call DrawFigure(ld,n,plane,lds,ns,square,1.0d0)
!
   call IntData2pgm(ld,n,plane,"plane")
!   
   stop
end program prot

subroutine IntData2pgm(s1,s2,idata,name)
!  Simple subroutine to dump integer data in a PGM format
   implicit none
   INTEGER, INTENT(IN) :: s1, s2
   INTEGER, DIMENSION(s1,s2), intent(in) :: idata
   character(*), intent(in) :: name
!
   integer :: i, j, ouni
   integer :: vp, vs
   real :: rmin, rmax
   character(80) :: fname
!
!  Write on unit 700 with  PGM format
   ouni = 700
   fname = TRIM(name)//".pgm"
   open(ouni,file=trim(fname),status="replace",iostat=vs)
   if ( vs /= 0 ) then
      print*,"!!!! Error write access to file ",trim(fname)
   end if
!  Magic code
   write(ouni,"(a)") "P2"
!  Dimensions
   write(ouni,*) s1, s2
!  Maximum value
   write(ouni,"(a)") "255"
!  Values from 0 to 255
   rmin = MINVAL(idata); rmax = MAXVAL(idata)
   vs = 0
   do i = 1, s1
      do j = 1, s2
         vp = INT ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) )
         vs = vs + 1
         write(ouni,"(i4)",advance="no") vp
         if (vs >= 10 ) then
            write(ouni,"(a)") " "
            vs = 0
         end if
      end do
      write(ouni,"(a)") " "
      vs = 0
   end do
   close(ouni)
   return
end subroutine IntData2pgm

subroutine DrawFigure(ld,n,plane,lds,ns,square,opaq)
   implicit none
   
   integer, intent(in) :: ld, n, lds, ns
   integer, intent(in out), dimension(ld,n) :: plane
   real(8), intent(in), dimension(lds,ns,2) :: square
   real(8), intent(in) :: opaq
   
   integer :: vp, i, j, ix, iy, ii, jj
   
   vp = int(256.0d0*opaq)
   if ( vp > 255 ) vp = 255
   
   do i = 1, lds
      do j = 1, ns
	      ix = square(i,j,1)+n/2; iy = square(i,j,2)+n/2
	      if ( ix > 0 .and. ix < ld ) then
	      	  if ( iy > 0 .and. iy < n ) then
	      	  	  ! draw a bold line
	      	  	  do ii = floor(square(i,j,1)+dble(n/2)), ceiling(square(i,j,1)+dble(n/2))
        	  	  	do jj = floor(square(i,j,2)+dble(n/2)), ceiling(square(i,j,2)+dble(n/2))
            	  	  plane(ii,jj) = vp
            	  	end do
            	  end do
	      	  end if
	      end if
      end do
   end do
   
   return
end subroutine DrawFigure