program rotation
!
! Rotate the external B field
!
       implicit double precision (a-h,o-z)
!
!      version = 0.01
!
       integer mxzgr, mxrgr
       character(80) title
       parameter ( mxxg=100, mxyg=100, mxzg=500 )
       common/XYSGRID2/xgr(mxxg),ygr(mxyg),sgr(mxzg)
       open(unit=13,file='simple_sol3D.dat',status='old',iostat=ioc)
       open(unit=23,file='rotated.dat',status='unknown',iostat=ioc1)
       if( ioc .ne. 0 ) then
       write(*,*) 'Cant open file ','simple_sol3D.dat'
        write(2,*) 'Cant open file ','simple_sol3D.dat'
        end if
       read(13,'(a80)',iostat=ioc) title
       write(2,'(a80)') title
       write(23,'(a80)',iostat=ioc1) title
       read(13,*,iostat=ioc) nxgr,nygr,nsgr
       write(23,*,iostat=ioc1) nxgr,nygr,nsgr
       if( ioc .ne. 0 ) then
        write(2,*)'Error reading nxgr,nygr,nsgr'
        stop
       end if
       if( ioc1 .ne. 0 ) then
        write(2,*)'Error writing nxgr,nygr,nsgr'
        stop
       end if
      if( nxgr.gt.mxxg .or. nygr.gt.mxyg .or. nsgr.gt.mxzg ) then
        write(2,*)'Dimensions exceed max limit for 3D field'
        stop
      end if
      read(13,*) (xgr(j),j=1,nxgr)
      read(13,*) (ygr(j),j=1,nygr)
      read(13,*) (sgr(j),j=1,nsgr)
      write(23,*) (xgr(j),j=1,nxgr)
      write(23,*) (ygr(j),j=1,nygr)
      write(23,*) (sgr(j),j=1,nsgr)
!
      nrecords = nxgr*nygr*nsgr
      dxgr = xgr(2) - xgr(1)
      dygr = ygr(2) - ygr(1)
      dsgr = sgr(2) - sgr(1)
       piover2=asin(1.d0)
       write(6,*)'enter: B scale'
       read(5,*)cscale
       write(6,*)'enter: rotation angle'
       read(5,*)fi
       fix=fi*piover2/90.0d0
       write(6,*)'enter: direction angle of axis of rotation n=(cos(theta),sin(theta),0)'
       read(5,*)theta
       THX=theta*piover2/90.d0
      do i=1,nrecords
      read(13,*,iostat=ioc)ix,iy,iz,bx,by,bz
       if( ioc .ne. 0 ) then
        write(2,*)'Error reading field data file'
        stop
       end if
       bxgr = bx * cscale
       bygr = by * cscale
       bsgr = bz * cscale
! rotate the field
       bx= bxgr*(cos(fix)+cos(THX)**2*(1-cos(fix)))+bygr*(sin(THX)*cos(THX)*(1-cos(fix)))+bsgr*sin(THX)*sin(fix)
       by = bxgr*(sin(THX)*cos(THX)*(1-cos(fix)))+bygr*(cos(fix)+sin(THX)**2*(1-cos(fix)))+bsgr*cos(THX)*sin(fix)
       bz = bxgr*sin(THX)*sin(fix)- bygr*cos(THX)*sin(fix)+bsgr*cos(fix)
        write(23,*,iostat=ioc1)ix,iy,iz,bx,by,bz
       if( ioc1 .ne. 0 ) then
        write(2,*)'Error writing field data file'
        stop
       end if

       end do
       
!
       close(13)
       close(23)
       stop
       end