      subroutine ellipse(xc,yc,major,minor,ang)
      real xc,yc,major,minor,ang
      parameter (n=41)
C      integer ia(idim,jdim)
C major and minor are like radii
C
C ang is in radians

      integer i,j
      real x(n),y(n),xtrans(n),ytrans(n),x0,y0,pi
      real xrot(n),yrot(n),qang
C      write(*,*) major,minor,ang
      pi=4.0*atan(1.0)
C define coordinate of rectangle
      x0=major
      y0=minor

C      write(*,'(a,2f7.4)') 'x0,y0 ',x0,y0
      do i=1,n
         x(i)=x0*cos(2.0*pi*float(i)/float(n)) !+xc
         y(i)=y0*sin(2.0*pi*float(i)/float(n)) !+yc
         xtrans(i)=x(i)*cos(ang)-y(i)*sin(ang)
         ytrans(i)=x(i)*sin(ang)+y(i)*cos(ang)
      enddo

C xtrans and ytrans are the ellipse boundaries in Qx space in AA^-1
C units
C
C the ellipse should be rotated by the angle of Q relative to Zeta

      qang=atan(yc/xc)

      do i=1,n
         xrot(i)=cos(qang)*xtrans(i)-sin(qang)*ytrans(i)
         yrot(i)=sin(qang)*xtrans(i)+cos(qang)*ytrans(i)
      enddo

      do i=1,n
         xrot(i)=xc+xrot(i)
         yrot(i)=yc+yrot(i)
      enddo
C      call PGPIXL (IA, IDIM, JDIM, 1, 1, idim, jdim, 
C     1                   X1, X2, Y1, Y2)
      call pgpoly(n,xrot,yrot)
      return
      end

C-----------------------------------------------------------------

      subroutine ellipseproj(xc,yc,arm,um)

      double precision arm(3,3),um(3,3)
      real xc,yc,major,minor,ang
      parameter (n=41)
C      integer ia(idim,jdim)
C major and minor are like radii
C
C ang is in radians

      integer i,j
      real yy(3,n)
      real x(n),y(n),xtrans(n),ytrans(n),x0,y0,pi
      real xrot(n),yrot(n),qang
C      write(*,*) major,minor,ang
      pi=4.0*atan(1.0)
C define coordinate of rectangle
      x0=sqrt(2.0*log(2.0)/arm(1,1))
      y0=sqrt(2.0*log(2.0)/arm(3,3))
C

C principle axes coordinates
      do i=1,n
         yy(1,i)=x0*cos(2.0*pi*float(i)/float(n)) !+xc
         yy(2,i)=0.0
         yy(3,i)=y0*sin(2.0*pi*float(i)/float(n)) !+yc
C         xtrans(i)=x(i)*cos(ang)-y(i)*sin(ang)
C         ytrans(i)=x(i)*sin(ang)+y(i)*cos(ang)
      enddo
C transform to q's

      do i=1,n
         xtrans(i)=yy(1,i)*um(1,1)+yy(3,i)*um(1,3)
         ytrans(i)=yy(1,i)*um(2,1)+yy(3,i)*um(2,3)
      enddo

C xtrans and ytrans are the ellipse boundaries in Qx space in AA^-1
C units
C
C the ellipse should be rotated by the angle of Q relative to Zeta

      qang=atan(yc/xc)

      do i=1,n
         xrot(i)=cos(qang)*xtrans(i)-sin(qang)*ytrans(i)
         yrot(i)=sin(qang)*xtrans(i)+cos(qang)*ytrans(i)
      enddo

      do i=1,n
         xrot(i)=xc+xrot(i)
         yrot(i)=yc+yrot(i)
      enddo
C      call PGPIXL (IA, IDIM, JDIM, 1, 1, idim, jdim, 
C     1                   X1, X2, Y1, Y2)
      call pgpoly(n,xrot,yrot)
      return
      end








