MODULE Prototypes
   IMPLICIT NONE
!
   INTERFACE
   !
   	 SUBROUTINE DROTG(DA,DB,C,S)
   	 !  DROTG construct givens plane rotation.
        IMPLICIT NONE
        REAL(8), INTENT(IN OUT) :: DA, DB
        REAL(8), INTENT(OUT) :: C, S
     END SUBROUTINE DROTG
     
     SUBROUTINE DROT(N, DX, INCX, DY, INCY, C, S)
     !  DROT applies a plane rotation.
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: N, INCX, INCY 
        REAL(8), INTENT(IN) :: C, S
        REAL(8), INTENT(IN OUT) :: DX, DY
     END SUBROUTINE DROT
     
   END INTERFACE
   
END MODULE Prototypes
   	   
 


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, alfa, alfa_rad, div
   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*,"Rotation angle in degrees? "
   read(*,*) alfa
   call xy_compute(da,db,alfa)
   
!
   !!! code here	! get parameters C, S
   !!! code here to realize rotation
   
!       
!	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


subroutine xy_compute(da,db,alfa)
!     Computes da, db given an angle in degrees
      implicit none
      real(8), intent(out) :: da, db
      real(8), intent(in out) :: alfa
      
      real(8) :: div, alfa_rad
      
      if ( alfa > 360.0d0 ) then
      	  div = int(alfa/360.0d0)
      	  alfa = alfa - div*360.0d0
      end if
      if ( abs(alfa) == 90.0d0 ) then
      	  da = 0.0d0; db = 1.0d0; alfa_rad = 2.0d0*atan(1.0d0)
      else if ( abs(alfa) == 270.0d0 ) then
      	  da = 0.0d0; db = -1.0d0; alfa_rad = 6.0d0*atan(1.0d0)
      else    	  
         alfa_rad = alfa*8.0d0*atan(1.0d0)/360.0d0
         da = 1.0; db = da*tan(alfa_rad)
      end if
      return
end subroutine xy_compute