ccc f77 ellipse_ha.f -o ellipse_ha.e -L/star/lib `cfitsio_link` `sla_link` -l pgplot -L/usr/lib64 -lX11 car_pltel.f -fno-second-underscore

cc       subroutine ellipse(ra_cen,dec_cen,x,y,rarray,darray,
cc     >magarray,major,minor,pa)
ccc    f77 ellipse.f -o ellipse.e -L/star/lib `fitsio_link` `sla_link`
ccccc `pgplot_link` car_pltel.f
       implicit none
       real*8 ra_cen,dec_cen,epoch,TEMPRA,TEMPDEC,DUMRA,DUMDEC,
     > sla_dsep,pi,sep,xi,eta,nulld,dra,ddec,gall,galb,dreslt
       real magarray,major,minor,pa,xsq(4),ysq(4),
     > x,y,plotx,ploty,axadd,ayadd,bxadd,byadd,xpt(40),ypt(40),
     > bright,faint,ell
       integer i,j,status,equinox,outequinox,purge,subset,npts,
     > blend,quality,form,wave,nstrt,jf,step,tempv,nrows,area,
     > naxis2,whichp,ifcor,total
       character*30 string
       character cstring*28,tstring*25,outepoch*17,comment*32
       character*5 pid,cobj
       character*14 device 
       character*1 dum,t
       character*6 numstring
       character*150 line
       character*28 title
       character*36 title2
      integer unit,readwrite,blocksize,hdutype,ntable,noobj
      integer felem,nelems,nullj,diameter,nfound,irow,colnum
      real nulle,density,radius
      character filename*40,nullstr*1,name*8,ttype(3)*10
      logical anynull
      double precision c1,c2
      real c3
      double precision c4,c5
      real c6,c7,c8,c9
      integer c10,c11,c12,c13
      real c14,c15,c16
      real c17,c18,c19,c20,c21
      integer length
      integer c22,c23,c24,c25,c26,c27,
     :c28,c29
      integer*2 c30,c32
      byte c31
      integer c33,c34,c35,c36,c37,
     :c38,c39,c40,c41,c42,c43,c44,c45,c46,c47

       pi=3.141592654d0
       nstrt=1
       read (5,33) pid
 33    format(a5)
       read (5,*) wave
       read (5,'(a)') string
       read (5,*) equinox
       read (5,*) x
       read (5,*) y
       read (5,*) outequinox
       read (5,'(a)') outepoch
       read (5,*) bright
       read (5,*) faint
       read (5,*) purge
       read (5,*) subset
       read (5,*) blend
       read (5,*) quality
       read (5,*) form
       read (5,*) whichp
       read (5,*) ifcor
       read (5,*) epoch
       read (5,*) noobj
       read (5,'(a)') device
        if (epoch.eq.0.0) epoch=1980.0

       radius=x
       if (y.eq.0.) then
       radius=y
       x=x*2.
       y=x
       end if
      

       call pgbegin(0,device,1,1)
       call pgslw(5)
      call pgsch (1.5)
       call pgenv(x*.6,-x*.6,-y*.6,y*.6,1,0)

c       call pgenv(-x/2,x/2.,-y/2.,y/2.,1,0)

         
       step=1

       if (equinox.eq.1 .or. equinox.eq. 2) then
       call sla_dafin (string,nstrt,ra_cen,jf)
       call sla_dafin (string,nstrt,dec_cen,jf)
       ra_cen=ra_cen*15.
       endif

       if (equinox.eq.0 .and. outequinox .eq. 1) then
       CALL sla_DFLTIN (STRING, NSTRT, gall, JF) 
       CALL sla_DFLTIN (STRING, NSTRT, galb, JF) 
       CALL sla_GE50 (gall*pi/180.d0, galB*pi/180.d0, ra_cen, dec_cen)
       end if

       if (equinox.eq.0 .and. outequinox .eq. 2) then
       CALL sla_DFLTIN (STRING, NSTRT, gall, JF) 
       CALL sla_DFLTIN (STRING, NSTRT, galb, JF) 
       CALL sla_galeq (gall*pi/180.d0, galB*pi/180.d0, ra_cen, dec_cen)
       end if
 
       if (outequinox.eq.1 .and. equinox.eq.2) then
       call sla_fk54z (ra_cen,dec_cen,epoch,tempra,tempdec,
     >dumra,dumdec)
       ra_cen=tempra
       dec_cen=tempdec
       end if
      if (outequinox.eq.2 .and. equinox.eq.1) then
       call sla_fk45z (ra_cen,dec_cen,epoch,tempra,tempdec)
       ra_cen=tempra
       dec_cen=tempdec
       end if

       call pgsfs(2)

       if (form.eq.1) then
       write (6,*)'<br>'
        if (noobj .gt. 2000) step=1+ nint(float(noobj)/2000.)
       if (step .le. 1) step=1
       open (unit=13, file=pid//'.ascii',status='old')
       do i=1,3
       read (13,10) dum
       end do
       do i=1,10000000
       
       nstrt=1

       read (13,22,end=99) cstring,
     :c3,c4,c5,c6,c7,c8,
     :c9,c10,c11,c12,c13,c14,c15,
     :c16,c17,c18,c19,c20,c21,c22,
     :c23,c24,c25,c26,c27,c28,c29,
     :c30,c31,c32
 22   format(a28,3x,f9.3,
     :2f12.6,4f9.3,4i11,4f9.4,4f10.3,4i10,4i12,i6,i4,i6)

       

       call sla_dafin (cstring,nstrt,dra,jf)
       call sla_dafin (cstring,nstrt,ddec,jf)
       

       
       dra=dra*15.
cccc       write (6,*)cstring,' ',dra,ddec,ra_cen,dec_cen
       pa=float(c30)
       pa=90.-pa
       pa=pa*pi/180.
       if (wave .eq. 1) then
       ell=c14
       area=c10
       end if
       if (wave .eq. 2) then 
       ell=c15
       area=c11
       end if

        major=sqrt(float(area)/(pi*(1.-ell)))
        minor=float(area)/(pi*major)
        major=major*.7/60.
        minor=minor*.7/60.

       
       call sla_DS2TP(dra,ddec,ra_cen,dec_cen,
     >xi,eta,j)
       plotx=xi*180.*60./pi
       ploty=eta*180.*60./pi



       call car_pltel(plotx,ploty,major,minor,pa,40,npts,
     >xpt,ypt,status)

       call pgpoly (npts,xpt,ypt)

cccc       write (6,*)' ',plotx,ploty
       do j=1,step-1
       read (13,10,end=99) dum
 10    format(a1)
       enddo   
       end do
       end if
       if (form.eq.2) then
       write (6,*)'<br>'
       open (unit=13, file=pid//'.tab',status='old')

       if (noobj .gt. 2000) step=1+ nint(float(noobj)/2000.)
       do i=1,15
       read (13,10) dum
       end do
       if (step .le. 1) step=1
       do i=1,10000000

       
       nstrt=1


       read (13,44,end=99) 
     :total,t,c1,t,c2,t,
     :c3,t,c4,t,c5,t,c6,t,c7,t,c8,t,
     :c9,t,c10,t,c11,t,c12,t,c13,t,c14,t,
     :c15,t,c16,t,c17,t,c18,t,c19,t,
     :c20,t,c21,t,c22,t,c23,t,c24,t,c25,t,
     :c26,t,c27,t,c28,t,c29,t,c30,t,c31,t,
     :c32
44    format(i6.6,a1,f11.7,a1,f11.7,a1,f9.3,a1,
     :2(f11.6,a1),4(f9.3,a1),4(i11,a1),4(f9.4,a1),4(f10.3,a1),4(i8,a1),
     :4(i12,a1),i6,a1,i4,a1,i6)

       dra=c1
       ddec=c2
       dra=dra*pi/180.
       ddec=ddec*pi/180.

       pa=float(c30)
       pa=90.-pa
       pa=pa*pi/180.
       if (wave .eq. 1) then
       ell=c14
       area=c10
       end if
       if (wave .eq. 2) then 
       ell=c15
       area=c11
       end if

        major=sqrt(float(area)/(pi*(1.-ell)))
        minor=float(area)/(pi*major)
        major=major*.7/60.
        minor=minor*.7/60.

       
       call sla_DS2TP(dra,ddec,ra_cen,dec_cen,
     >xi,eta,j)
       plotx=xi*180.*60./pi
       ploty=eta*180.*60./pi



       call car_pltel(plotx,ploty,major,minor,pa,40,npts,
     >xpt,ypt,status)

       call pgpoly (npts,xpt,ypt)
       write (6,*)' '
       do j=1,step-1
       read (13,10,end=99) dum
       enddo   
       end do
       end if

       if (form.eq.3) then
       write (6,*)'<br>'
        status=0

C     Get an unused Logical Unit Number to use to open the FITS file
c      call ftgiou(unit,status)
      unit=2
C     open the FITS file previously created by WRITEIMAGE
      filename=pid//'.fits'
      readwrite=0
      call ftopen(unit,filename,readwrite,blocksize,status)

      ntable=1


 5        call ftmrhd(unit,1,hdutype,status)
        



 6        call ftgkns(unit,'TTYPE',1,3,ttype,nfound,status)
ccccc          write(*,2000)ttype
2000      format(8x,3a10)

C         read the data, one row at a time, and print them out
          felem=1
          nelems=1
          nullstr=' '
          nullj=0
          nulle=0.
          nulld=0.
       call ftgkyj(unit,'NAXIS2',nrows,comment,status)
         
          noobj=nrows
        
        if (noobj .gt. 2000) step=1+nint(float(noobj)/2000.)
          do irow=1,nrows,step


 7            call ftgcvd(unit,1,irow,felem,nelems,nulld,
     >dra,anynull,status)
             dra=dra*pi/180.0d0

              
 8            call ftgcvd(unit,2,irow,felem,nelems,nulld,
     >ddec,anynull,status)
        ddec=ddec*pi/180.0d0

 9            call ftgcvi(unit,30,irow,felem,nelems,nullj,
     >c30,anynull,status)


       if (wave .eq. 1) then
           call ftgcvj(unit,10,irow,felem,nelems,nullj,
     >area,anynull,status)
           call ftgcve(unit,14,irow,felem,nelems,nullj,
     >ell,anynull,status)

       end if
       if (wave .eq. 2) then 
           call ftgcvj(unit,11,irow,felem,nelems,nullj,
     >area,anynull,status)
           call ftgcve(unit,15,irow,felem,nelems,nullj,
     >ell,anynull,status)

       end if


       pa=float(c30)


       pa=90.-pa
       pa=pa*pi/180.
            

        major=sqrt(float(area)/(pi*(1.-ell)))
        minor=float(area)/(pi*major)
        major=major*.7/60.
        minor=minor*.7/60.
       
       call sla_DS2TP(dra,ddec,ra_cen,dec_cen,
     >xi,eta,j)
       plotx=xi*180.*60./pi
       ploty=eta*180.*60./pi



       call car_pltel(plotx,ploty,major,minor,pa,40,npts,
     >xpt,ypt,status)

       call pgpoly (npts,xpt,ypt)
       write (6,*)' '

      
          end do
      

C     close the file and free the unit number
 11       call ftclos(unit, status)
      call ftfiou(unit, status)
    
C     check for any error, and if so print out error messages

       

       end if

 99    continue
      
     
       if (step .gt. 1) then
         write (6,*) ' ',noobj, ' objects were extracted<p>'
       write (6,*) ' For clarity and efficiency only 1 object in '
       write (6,*) step
       write (6,*) ' has been plotted in the chart below<p>'
       write (cobj,123) step
       title='1 in '//cobj//' objects plotted'
       else
       write (cobj,123) noobj
       title=cobj//' objects plotted'
       endif
 123   format(i5)
       call pglab('<- RA (arcmin)','Dec (arcmin) ->',title)
       call pgsls(4)
       if (noobj.lt.100) then
       if (radius.eq.0.0) then

       call pgcirc(0.,0.,x/2)
       else
       xsq(1)=-x/2
       ysq(1)=-y/2
       xsq(2)=-x/2
       ysq(2)=y/2
       xsq(3)=x/2
       ysq(3)=y/2
       xsq(4)=x/2
       ysq(4)=-y/2
       call pgpoly (4,xsq,ysq)
       endif
       end if
       title2=string
       call pgptxt(0.0,y*.65,0.0,0.5,title2)
       call pgsch(1.0)
       call pgptxt(0.0,y*.55,0.0,0.5,'N')
       call pgptxt(x*0.55,0.0,0.0,0.5,'E')

       call pgend
 
       end



      subroutine printerror(status)

C     Print out the FITSIO error messages to the user

      integer status
      character errtext*30,errmessage*80

C     check if status is OK (no error); if so, simply return
      if (status .le. 0)return

C     get the text string which describes the error
 1    call ftgerr(status,errtext)
      print *,'FITSIO Error Status =',status,': ',errtext

C     read and print out all the error messages on the FITSIO stack
 2    call ftgmsg(errmessage)
      do while (errmessage .ne. ' ')
          print *,errmessage
          call ftgmsg(errmessage)
      end do
      end





























