      subroutine ellipse(xc,yc,major,minor,ang)
      implicit none
c...	constants
      integer n
      parameter (n=41)
c...	passed variables
      real xc,yc,major,minor,ang

c...	local variables
      
C      integer ia(idim,jdim)
C major and minor are like radii
C
C ang is in radians

      integer i
c0111      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,rm,arm,um)
      implicit none
      include 'spur_color.inc'

c..	constants
      integer n
      parameter (n=41)

c...	passed variables
      double precision rm(4,4),arm(3,3),um(3,3)
      real xc,yc

c...	local variables      
c0111Unused      real major,minor,ang
      double precision ap,b,c,v2,v,p1,p2
c0111Unused     double precision alfa1,alfa2
      double precision fwhmt1
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,n)
      real xtrans(n,n),ytrans(n,n)
      real x0,y0, z0
      real pi
      real theta, phi
c0111Unused      real x(n),y(n)
      real etrans(n,n)
c0111Unused      real xrot(n),yrot(n),qang
      real temp1, temp2

C      write(*,*) major,minor,ang
      pi=4.0*atan(1.0)

      fwhmt1=1.0d0/dsqrt(8.0d0*dlog(2.0d0))
      AP=RM(1,1)-RM(1,3)*RM(1,3)/RM(3,3)
      B=RM(2,2)-RM(2,3)*RM(2,3)/RM(3,3)
      C=RM(1,3)*RM(2,3)/RM(3,3)-RM(1,2)
      V2=DATAN(2*C/(B-AP))
      V=V2/2
      P1=1/DSQRT(AP*DCOS(V)**2+B*DSIN(V)**2-C*DSIN(V2))/FWHMT1
      P2=1/DSQRT(AP*DSIN(V)**2+B*DCOS(V)**2+C*DSIN(V2))/FWHMT1

C define major axes

      x0=sqrt(2.0*log(2.0)/arm(1,1))
      y0=sqrt(2.0*log(2.0)/arm(2,2))
      z0=sqrt(2.0*log(2.0)/arm(3,3))

C parametrization of principle axes coordinates

      do i=1,n
         do j=1,n
            theta=2.0*pi*float(i)/float(n)
            phi=pi*float(j)/float(n)
            yy(1,i,j)=x0*cos(theta)*sin(phi)
            yy(2,i,j)=y0*sin(theta)*sin(phi)
            yy(3,i,j)=z0*cos(phi)
         enddo
      enddo

C transform to q's

      do i=1,n
         do j=1,n
            xtrans(i,j)=yy(1,i,j)*um(1,1)+yy(2,i,j)*um(1,2)+
     &           yy(3,i,j)*um(1,3)
            ytrans(i,j)=yy(1,i,j)*um(2,1)+yy(2,i,j)*um(2,2)+
     &           yy(3,i,j)*um(2,3)
            etrans(i,j)=yy(1,i,j)*um(3,1)+yy(2,i,j)*um(3,2)+
     &           yy(3,i,j)*um(3,3)
	    temp1 = xc+xtrans(i,j)
	    temp2 = yc+ytrans(i,j)
            call pgcirc(temp1, temp2, .002)
c01Nov            call pgcirc(xc+xtrans(i,j),yc+ytrans(i,j),.002)
         enddo
      enddo
C plot extremum points
      call pgsci(color_blue)
c01Dec      call pgsci(4)
      call pgsfs(2)
      temp1 = sngl(p1)
      temp2 = sngl(p2)
      call ellipse(xc,yc,temp1, temp2,sngl(v))
c0111      call ellipse(xc,yc,float(p1),float(p2),v)

C 1,0,0
      write(*,*) 'extrema ',xc+x0*um(1,1),yc+x0*um(1,2)
      temp1 = sngl(xc+x0*um(1,1))
      temp2 = sngl(yc+x0*um(1,2))
      call pgcirc( temp1, temp2, .1)
c0111      call pgcirc(xc+x0*um(1,1),yc+x0*um(1,2),.1)
C -1,0,0
      temp1 = sngl(xc-x0*um(1,1))
      temp2 = sngl(yc-x0*um(1,2))
      call pgcirc(temp1, temp2, .1)
c0111      call pgcirc(xc-x0*um(1,1),yc-x0*um(1,2),.1)


      return
C      call PGPIXL (IA, IDIM, JDIM, 1, 1, idim, jdim, 
C     1                   X1, X2, Y1, Y2)

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

c0111	Due to return, this line can not be reached
c0111      qang=atan(yc/xc)

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

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








