*+  CAR_PLTEL - Calculates array of points describing an ellipse
      SUBROUTINE car_pltel (X, Y, MAJAX, MINAX, THETA, SIZE, NPTS, 
     :             PX, PY, STATUS)
*    Description :
*     This subroutine calculates an array of points describing an 
*     ellipse centered on x,y with given semi-major and semi-minor
*     axes and orientation.
*    Invocation :
*      CALL CAR_PLTEL (X, Y, MAJAX, MINAX, THETA, SIZE, NPTS, PX, PY, 
*    :        STATUS)
*
*    Parameters :
*     X             = REAL (READ)
*            x value
*     Y             = REAL (READ)
*            y value
*     MAJAX         = REAL (READ)
*            Semi-major axis value
*     MINAX         = REAL (READ)
*            Semi-minor axis value
*     THETA         = REAL (READ)
*            Orientation value
*     SIZE          = INTEGER (READ)
*            maximum size of array describing ellipse
*     NPTS          = INTEGER (WRITE)
*            Size of array describing ellipse
*     PX(SIZE)      = INTEGER (WRITE)
*            Array of x values for plotting
*     PY(SIZE)      = INTEGER (WRITE)
*            Array of y values for plotting
*     STATUS        = INTEGER (UPATE)
*            Running status
*    Method :
*      If (major axis ~ 0 ) then
*        draw a dot
*      else if (minor axis ~ 0 ) then
*        draw a line with length equal to the major 
*         axis and with given orientation
*      else
*        for array of points,
*          compute an angle phi
*          compute radius of the ellipse at this angle
*          compute corresponding X & Y positions
*        end do
*        for array of points
*          rotate to desired inclination
*        end do
*      end if
*      end
*    Deficiencies :
*    Bugs :
*     None known.
*    Authors :
*     S K Leggett.  (ROE::SKL)
*    History :
*     8/2/88:  Original version.                        (ROE::SKL)
*    Type Definitions :
      IMPLICIT NONE
*    Global constants :
c      INCLUDE 'SAE_PAR'
c      INCLUDE 'adc_$err'
*    Import :
      INTEGER
     :  SIZE           ! maximum size of array describing ellipse
      REAL
     :  X,             ! x value
     :  Y,             ! y value
     :  MAJAX,         ! Semi-major axis value
     :  MINAX,         ! Semi-minor axis value
     :  THETA          ! Orientation value
*    Import-Export :
*    Export :
      INTEGER
     :  NPTS           ! Size of array of values describing ellipse
      REAL
     :  PX(SIZE),      ! Array of x values for plotting
     :  PY(SIZE)       ! Array of y values for plotting
*    Status :
      INTEGER
     :  STATUS         ! Running status.
*    External references :
*    Global variables :
*     <any INCLUDE files for global variables held in named COMMON>
*    Local Constants :
      REAL PI
*      PARAMETER (PI=3.1415927)
*    Local variables :
      REAL 
     :  PHI,           ! Angle for calculating ellipse positions
     :  DEPHI,         ! Increment on angle 
     :  R,             ! Radius for calculating ellipse positions
     :  TEMPX,         ! Temporary value of PX
     :  TEMPY          ! Temporary value of PY
      INTEGER 
     :  INDEX          ! Counter
*    Internal References :
*    Local data :
*-

      pi=3.1415927
D     PRINT3000, STATUS
D3000 FORMAT( / 1X, 'CAR_PLTEL on entry;' /
D    :  2X, 'STATUS : ', I7 /)



c      IF (STATUS .EQ. SAI__OK) THEN




          IF ( (2.0*MAJAX) .LE. 1.0E-10) THEN


*           The major axis is sensibly zero. Draw a dot.

              NPTS  = 1
              PX(1) = X
              PY(1) = Y


          ELSE IF ( (2.0*MINAX) .LE. 1.0E-10) THEN


*           The eccentricity is sensibly 1.0. Draw a line.


              NPTS  = 2
              PX(1) = X - ( MAJAX * COS(THETA) )
              PY(1) = Y - ( MAJAX * SIN(THETA) )
              PX(2) = X + ( MAJAX * COS(THETA) )
              PY(2) = Y + ( MAJAX * SIN(THETA) )




          ELSE


*           Now do the case where a real ellipse has been input.



              NPTS = SIZE
   

*           Compute points along a horizontal, centred at origin

              DEPHI = (2.0*PI) / REAL( (SIZE-1) )

              DO INDEX = 1, SIZE

                  PHI = DEPHI * REAL(INDEX-1)

                  R = SQRT( 
     :              ((MAJAX*MINAX)**2) /  
     :              ( ((MAJAX*SIN(PHI))**2) + ((MINAX*COS(PHI))**2) ) )

                  PX(INDEX) = R * COS(PHI)
                  PY(INDEX) = R * SIN(PHI)

              END DO


*           Rotate through an angle theta, & translate to x,y

              DO INDEX = 1, SIZE

                  TEMPX = PX(INDEX)
                  TEMPY = PY(INDEX)

                  PX(INDEX) = X + ( TEMPX * COS(THETA) ) 
     :              - ( TEMPY * SIN(THETA) )
                  PY(INDEX) = Y + ( TEMPX * SIN(THETA) )
     :              + ( TEMPY * COS(THETA) )

              END DO



          END IF



c      END IF

 
D     PRINT4000, STATUS
D4000 FORMAT( / 1X, 'CAR_PLTEL on exit;' /
D    :  2X, 'STATUS : ', I7 /)


      END
















