C*****************************************
      program maskgen4
c
c	INPUT FORMAT:

c	col(1) Object Number: integer; any order; can be output from
c	       sextract or focas; include alignment stars.

c	col(2), (3): Offset from the nominal mask center in arcsec
c	             e/w and n/s respectively: e is + and n is +


ccccccccccc     REBECCA NOTES THAT PRECESSION WILL CAUSE THE APPARENT PA
c		TO SHIFT, SO SHOULD CALCULATE THIS SHIFT, GIVEN THE EPOCH
c		OF THE POSITIONS


c	col 4 :  priority number(floating) > 0; LOWEST #s have HIGHEST priority
c                the priority for alignment stars should be given a
c                priority of -1.0.  You can put magnitude in this column (except
c                for alignment stars), to prioritize the objects by brightness.
c
c	col 5: magnitude 
c
c       col 6: Position Angle of slit, for tilted slit (optional).  Note that
c              180>PA>0.  If E-W slits are used (recommended) then CASS RING angle 
c              is 270, and the most effective position angles are 135>PA>45.
c
c              Note that tilted slits have an increased width in the dispersion
c              direction, resulting in a loss of resolution, relative to untilted slits.
c
                     

c the program will output: ...

c		1) A plot on screen of selected masks
c
c   optional    2) A postscript file named "photomask.ps" necessary for producing a
c                    photographic slit mask.  Numbers can be included to make an 11X17
c                    identification chart (ioption 6).
c
c               3) A postscript file named "pgplot.ps" for identification of the
c                    slit mask.
c
c		4) Output file containing the alignment stars and all the other objects
c                    NOT in the present mask file.
c
c   optional    5) An output file containg all object used in the current mask (ioption 5).
c 
c               6) Two output files, username.ps1 and username.ps2, which contain the 
c                    GM code instructions for the milling machine to cut a metal slit mask.
c


c
c It is necessary to have a local file named "maskgen.common.blocks" when compiling.
c The file contains the following statements:
c
c--------------start of maskgen.common.blocks---------------------------------
c
c      integer kmax
c      parameter (kmax=4000)  ! maximum number of objects in a mask list
c
c
c      common/iindata/ inumin,inumobj
c      integer inumin,inumobj(kmax)
c
c      common/islit/ iobjproj,iflag,kcount,icount,iscale
c      integer iobjproj(kmax),iflag(kmax),kcount(kmax),icount,iscale
c
c      common/iparams/ isense,ixsign,iysign,iwarn,ioffset,kpri,iwfccd,
c    %                iskip
c      integer isense,ixsign,iysign,iwarn,ioffset,kpri,iwfccd,iskip
c
c      common/rindata/xobj,yobj,priority,xmag,pa,wslit,slitposw,slitpose
c      real xobj(kmax),yobj(kmax),priority(kmax),xmag(kmax),pa(kmax)
c      real wslit(kmax),slitposw(kmax),slitpose(kmax)
c
c      common/rslit/ xproj,yproj,xstart,xend,ystart,yend,slitposwproj,
c    %              slitposeproj,x0,y0,x0unproj,y0unproj,wproj,lproj,
c    %              paproj,least,lwest,percent,projpriority
c      real xproj(kmax),yproj(kmax),xstart(kmax),xend(kmax),ystart(kmax)
c      real yend(kmax),slitposwproj(kmax),slitposeproj(kmax),x0,y0
c      real x0unproj,y0unproj,wproj(kmax),lproj(kmax),least(kmax)
c      real paproj(kmax),lwest(kmax),percent(kmax),projpriority(kmax)
c
c      common/rparams/ platescale,scale,starwidth,xsize,ysize,ystack,
c    %                endpix,endslit,slitmin,enlarge,xoff,yoff,globslit,
c    %                thinrot,theta,thetad,cosmic_scale
c      real platescale,scale,starwidth,xsize,ysize,ystack,endpix,endslit
c      real slitmin,enlarge,xoff,yoff,globslit,thinrot,theta,thetad,cosmic_scale
c
c      common/lparams/ auto_priority, tilted_slits, holes, stacked, cosmic
c      logical auto_priority, tilted_slits, holes, stacked, cosmic
c
c      common/cparams/ title,infile
c      character*32 title,infile*64
c
c      common/rcoords/ rah,ram,ras,decd,decm,decs,epoch,rahproj,ramproj,
c    %                 rasproj,decdproj,decmproj,decsproj
c      real rah,ram,ras,decd,decm,decs,epoch,rahproj,ramproj,rasproj
c      real decdproj,decmproj,decsproj
c
c      common/lcoords/coords
c      logical coords
c
c--------------end of maskgen.common.blocks-----------------------------------
c


      include 'maskgen.common.blocks'

      data scale,cosmic_scale/10.91585,11.81/        
      data starwidth/8.0/  ! width of the alignment star boxes in arcsec
      data xsize,ysize,ystack/1500.0,1080.0,360.0/ !the limits on field in arcsec
      data endpix,endslit,slitmin,globslit/4.0,2.0,10.0,1.5/

      data enlarge/3.75/  ! enlargement factor
      data xoff,yoff/0.0,0.0/ ! intitialize some parameters...
      data holes,auto_priority,tilted_slits,coords/.false.,.true.,.false.,.false./
      data stacked,cosmic/.false.,.false./
      data isense,iscale,iwfccd,iskip,ioffset,kpri/0,0,0,0,0,0/
        
 

      call query_user         ! enter input file name, title and set parameters 

      call read_data

      call select_all_objects(iflag,inumin)   !initialize all objects to be considered

      call create_slit_mask                   ! select objects and create slit mask information

      call make_xwindow_plot(xsize,ysize,scale)  ! plot current slits on screen with PGPLOT

      call interact_with_user

      call produce_output_files

      call pgend
      stop  ' NORMAL STOP '
      end





      subroutine interact_with_user
c
c allow user to modify parameters and select options interactively
c
      include 'maskgen.common.blocks'

      logical numbers
      integer ioption

      numbers = .false.

      do while(.true.)

         call get_commands(ioption)

         if (ioption .eq. 1) then
            call change_cass_ring_pa
            call select_all_objects(iflag,inumin)
         elseif (ioption .eq. 2) then 
            call delete_obj
         elseif (ioption .eq. 3) then 
            call insert_obj
         elseif (ioption .eq. 4) then 
            continue                  ! replot current parameters/objects on screen
         elseif (ioption .eq. 5) then 
            call print_current_list
         elseif (ioption .eq. 6) then               ! make the photomask.ps file
            call prompt_for_numbers(numbers)
            dxpaper = 16.5
            dypaper = 10.5
            call make_postscript_file(dxpaper,dypaper,numbers)
            call rename_photomask_pgplot_file
         elseif (ioption .eq. 7) then 
            call move_slit_boundaries(inumobj,inumin,slitposw,slitpose)
         elseif (ioption .eq. 8) then 
            call select_all_objects(iflag,inumin)
         elseif (ioption .eq. 9) then 
            call exclude_all_objects(iflag,inumin)
         elseif (ioption .eq. 10) then 
            call change_offsets
            call select_all_objects(iflag,inumin)
         elseif (ioption .eq. 11) then 
            return
         else
            write(6, '(a)')'Enter value 1-11 only:'
            call get_commands(ioption)
         endif
         if (ioption.le.6 .or. ioption.eq.10) then
	    call create_slit_mask
            call make_xwindow_plot(xsize,ysize,scale)
         endif
      enddo

      return
      end







      subroutine multisort(n,arr,b1rr,b2rr,b3rr,nb4rr,b5rr,b6rr,nb7rr,
     %                     b8rr)
      parameter (kmax=4000)  ! maximum number of objects in a mask list
c
c
c               sort the objects into ascending x order, while keeping the other arrays
c               in the associated positions
c               subroutine is *adapted* from numerical recipes;
c               there are now 7 dependent vectors (yproj,wproj,iobjproj,
c               slitposeproj,slitposwproj,KCOUNT,projpriority) on the one sorting
c               vector (xproj); when the ascii string is added, this will be sorted
c               using a separate CHARSORT routine; a "spare" array for xproj must be
c               created for this.
c
c               ... NOTE THAT projpriority has been changed from fixed to
c               ... floating point: March 1, 1998
c

      real arr(kmax),b1rr(kmax),b2rr(kmax),b3rr(kmax)
      real b5rr(kmax),b6rr(kmax),b8rr(kmax)
      integer nb4rr(kmax),nb7rr(kmax)

      do j=2,n
         a=arr(j)
         b1=b1rr(j)
         b2=b2rr(j)
         b3=b3rr(j)
         nb4=nb4rr(j)
         b5=b5rr(j)
         b6=b6rr(j)
	 nb7=nb7rr(j)
         b8=b8rr(j)
         do i=j-1,1,-1
            if (arr(i) .le. a) goto 10
            arr(i+1)=arr(i)
            b1rr(i+1)=b1rr(i)
            b2rr(i+1)=b2rr(i)
            b3rr(i+1)=b3rr(i)
            nb4rr(i+1)=nb4rr(i)
            b5rr(i+1)=b5rr(i)
            b6rr(i+1)=b6rr(i)
            nb7rr(i+1)=nb7rr(i)
            b8rr(i+1)=b8rr(i)
          enddo
         i=0
 10      arr(i+1)=a
         b1rr(i+1)=b1
         b2rr(i+1)=b2
         b3rr(i+1)=b3
         nb4rr(i+1)=nb4
         b5rr(i+1)=b5
         b6rr(i+1)=b6
         nb7rr(i+1)=nb7
         b8rr(i+1)=b8
      enddo
      return
      end
C****************************** 
c       SUBROUTINE MVPG   !...this doesnt work for the numbered version of the ps file
c		... read in a mask file pg plot, insert the special
c		... line for 11 x 17 printing, and rename the file
c		... to your choice.
c	logical logic
c	character buf*1000,infilename*60, ext*3,outfilename*63
c	character bufout*1000
c	character INSERT*16,LAST*5,EXTRALINE*35
c	parameter(ext = '.ps')
c	parameter(INSERT ='%%BeginPageSetup')
c	parameter(LAST= '%%EOF')
c	parameter(EXTRALINE= 'statusdict begin 0 setpapertray end')
c
c	
c	open(unit=10,file='pgplot.ps',status='old')
c	print *, ' enter name of output .ps file(dont include .ps)'
c	read '(a)', infilename
c	mm=lnblnk(infilename)
c	mm1=mm+1
c	mm3=mm+3
c	outfilename(1:mm) = infilename
c	outfilename(mm1:mm3) = ext
c	print '(a)', outfilename
c	open(unit=20,file=outfilename,status='new')
c	logic = .true.
c	do while(logic)
c	    read(10,'(a)') buf
c	    mm=lnblnk(buf)
c	    bufout=buf(1:mm)
c	    write(20,'(a)') buf(1:mm)
c	    if(buf.eq.INSERT) logic = .false.
c	enddo
c	write(20,'(a)') EXTRALINE
c	logic = .true.
c	do while(logic)
c	    read(10,'(a)') buf
c	    mm=lnblnk(buf)
c	    bufout=buf(1:mm)
c	    write(20,'(a)') buf(1:mm)
c	    if(buf.eq.LAST) logic = .false.
c	enddo
c	RETURN
c	end
C****************************************************
 
        subroutine insert_obj
c 
c AMcW 16/11/98
c
c Insert a slit/object; but first check for conflicts and
c suggest deletions if necessary
c

      include 'maskgen.common.blocks'

      parameter (wfpix = 0.77402)  ! the wfccd pixel size in arcsec

      integer found_object

c
c read object id # and find place in the list
c
      iadd = 0
      do while (iadd.eq.0)
         write(6,'(35h Enter number of object to include:)')
         read(5,*)iinc
         iadd = found_object(iinc)
         if (iadd.eq.0) then
            write(6,'(45h No object found with that number, try again:)')
         endif
      enddo
c
c include the new object
c
      iflag(iadd) = 0
      icount = icount + 1
      xproj(icount)=xobj(iadd)*cos(theta) + yobj(iadd)*sin(theta)
c
c find the adjacent slits
c
      xlow = -0.5*xsize
      xup  =  0.5*xsize
      do j = 1,icount-1
         xdist =  xproj(j) - xproj(icount)
         if (xdist.lt.0.0) then
            if (xdist.ge.xlow) then
               ilow = kcount(j)
               xlow = xdist
            endif
         else
            if (xdist.le.xup) then
               iup = kcount(j)
               xup = xdist
            endif
         endif
      enddo
c
c now find conflicts with adjacent slits and suggest deletions
c

c stars too close
      dmin = endpix*wfpix + 2.0*endslit
      slength = xup - xlow - dmin
      if (-xlow.lt.dmin) then
         write(6,'(22h Conflict with object ,i10)')inumobj(ilow)
         call delete(inumobj(ilow),iflag(ilow))
      endif
      if (xup.lt.dmin) then
         write(6,'(22h Conflict with object ,i10)')inumobj(iup)
         call delete(inumobj(iup),iflag(iup))
      endif

c slit too small
      if (iflag(ilow).ne.0 .and. iflag(iup).ne.0 .and. slength.lt.slitmin) then
         write(6,'(51h Undersize slit: delete an adjacent slit suggested )')
         write(6,'(24h Adjacent slit ids are: ,2i10)')inumobj(ilow),inumobj(iup)
         call delete(inumobj(ilow),iflag(ilow))
         call delete(inumobj(iup),iflag(iup))
      endif
      return
      end



      integer function found_object(iinc)

      include 'maskgen.common.blocks'

      found_object = 0
      do i=1,inumin
         if (inumobj(i) .eq. iinc) then
            found_object = i
            return
         endif
      enddo
      return
      end



      subroutine delete_obj

      include 'maskgen.common.blocks'

      integer found_object,iobj,idel

      iobj = 0
      do while (.true.)
         write(6,'(34h Enter number of object to delete:)')
         read(5,*)idel
         iobj = found_object(idel)
         if (iobj.eq.0) then
            write(6,'(45h No object found with that number, try again:)')
         else
c           iflag(idel) = 1
            iflag(iobj) = 1
            return
         endif
      enddo
      return
      end





      subroutine delete(inum,iflag)

      character*1 ans

      write(6,'(30h Do you wish to delete object ,i10)')inum
      read(5,'(a1)')ans
      if (ans.eq.'y' .or. ans.eq.'Y') then
         iflag = 1
         write(6,'(16h Deleted object ,i10)')inum
      endif
      return
      end








        subroutine query_user

        include 'maskgen.common.blocks'

        call get_input_file_name(infile)
        call get_title(title)
        call get_coords
        call set_alignment_warning
        call change_defaults
        call change_cass_ring_pa
        call set_global_slitwidth

        return
        end




      subroutine read_data

      include 'maskgen.common.blocks'

      character*64 infile

      open(unit=11,file=infile,form='formatted',status='old')
      inumin=0

      if (iskip .gt. 0) then
         do i=1,iskip
            read(11,*,end=1001,err=1001)
         enddo
      endif

      do I=1,kmax   
         if (tilted_slits) then
            read(11,*,end=1001,err=1001) inumobj(i),xobj(i),
     %           yobj(i),priority(i),xmag(i),pa(i)
         else
            read(11,*,end=1001,err=1001) inumobj(i),xobj(i),
     %           yobj(i),priority(i),xmag(i)
         endif
         wslit(i) = globslit
         inumin=inumin+1
         if(inumin.ge.kmax) stop ' too many entries in input list'

         if (isense .eq. 1) then
                if (ixsign .eq. 1) then
                    xobj(i) = -xobj(i)
                endif
                if (iysign .eq. 1) then
                    yobj(i) = -yobj(i)
                endif
         endif
         if (iscale .eq. 1) then
            xobj(i)=xobj(i)*platescale
            yobj(i)=yobj(i)*platescale
         endif
         slitpose(i)=0.
         slitposw(i)=0.
c		... if this object is a star, set its slitwidth to "starwidth'
         if(priority(i).lt.0.0) wslit(i) = starwidth
      enddo 
 1001 close(unit=11,status='keep')

      return
      end




      subroutine make_new_xylist
c
c compute projected x,y positions of un-deleted objects and centroid
c

      include 'maskgen.common.blocks'

      icount=0
      x0=0.   
      y0=0.      ! intialize the centroid calculations... but normally we do not CHANGE
c                    ... the offset, but calculate the centroid: this is in the ROTATED
c                    ... system
c		     ... NOTE: x always increase from left to right on the pgplot and
c		     ...       y always increases from top to bottom; 
      



      do i=1,inumin     !  ... from the full input list, keep only the objects
c			      ... with iflag(i) = 0; icount is the running index for
c	                      ... the working set; rotate for each such object	
c		     IFLAG(IFULL) =1  means an objects has been dropped from the active
c			(working) list; =0 means it is still being considered for inclusion
c			in the final slit list...
	      

        if (iflag(i) .eq. 0) then
            icount=icount+1
            kcount(icount) = i   ! kcount(m) is the full count # corresponding
c				  ... to the mth object in the working vector list
            xproj(icount)    = xobj(i)*cos(theta) + yobj(i)*sin(theta)
            yproj(icount)    = yobj(i)*cos(theta) - xobj(i)*sin(theta)
            wproj(icount)    = wslit(i)
            paproj(icount)   = pa(i)
            iobjproj(icount) = inumobj(i)
            slitposeproj(icount)=slitpose(i)
            slitposwproj(icount)=slitposw(i)
            projpriority(icount) = priority(i)
            x0 = x0 + xproj(icount)
            y0 = y0 + yproj(icount)

        endif     ! endif for storing the unsorted, working vectors         
      enddo     !   end do for read in of working vectors 



c
c                       ...centroid of current working object list is (x0,y0)
c                       ...these centroids are in the current project of the (x,y) coordinates
      x0 = x0 / real(icount)
      y0 = y0 / real(icount)

c
c               ...the desired (x0,y0) are relative to the input (x=0,y=0) *non-projected*
c               ...coordinates, so derotate the coordinates of the centroid
c

      x0unproj=x0*cos(theta) - y0*sin(theta)
      y0unproj=y0*cos(theta) + x0*sin(theta)

      write(6,'(a)')'x0,y0'
      write(6,*)x0,y0
      write(6,'(a)')'x0unproj,y0unproj,theta'
      write(6,*) x0unproj,y0unproj,theta

      return
      end




      subroutine sort_slits
c
c sort slit information by ascending xproj
c
      include 'maskgen.common.blocks'

      call multisort(icount,xproj,yproj,paproj,wproj,iobjproj,slitposeproj,
     %                slitposwproj,kcount,projpriority)
      return
      end





        subroutine shift_xy

        include 'maskgen.common.blocks'

        integer i

c
c 		   ... transform all (x,y) positions relative to the centroid
c		   ...--normally bypass this
c
c                 .....option to bypass any automatic offsets


c	... if ioffset = 1, xoff,yoff only, =2, x0,y0 only/=3 xoff,yoff,x0,y0

        shiftx = 0.0
        shifty = 0.0

	if(ioffset.gt.0) then    !   ...make an offset correction:
	   if(ioffset.eq.1) then
              shiftx = -xoff
              shifty = -yoff
	   else if(ioffset.eq.2) then
              shiftx = -x0 
              shifty = -y0 
	   else if(ioffset.eq.3) then
              shiftx = -xoff - x0
              shifty = -yoff - y0
	   endif

           do i=1,icount
              xproj(i) = xproj(i)  + shiftx
              yproj(i) = yproj(i)  + shifty
           enddo
        endif   ! endif for auto-offset bypass

        call compute_shifted_coords(shiftx,shifty)

        return
        end




      subroutine compute_shifted_coords(shiftx,shifty)

      parameter (raddeg= 57.29577951)
      include 'maskgen.common.blocks'

      real shiftx, shifty

      if (.not.coords) return

c compute unprojected ra and dec shifts in arc seconds

      ra_shift  = shiftx*cos(theta) - shifty*sin(theta)
      dec_shift = shifty*cos(theta) + shiftx*sin(theta)

c convert to degrees

      call hms_to_deg(rah,ram,ras,ra_deg)
      call dms_to_deg(decd,decm,decs,dec_deg)

c compute new ra and dec

      dec_proj = dec_deg - dec_shift/3600.
      ra_proj  = ra_deg  - ra_shift/(3600.*cos(dec_deg/raddeg))

c convert to hms and dms system

      call deg_to_hms(ra_proj,rahproj,ramproj,rasproj)
      call deg_to_dms(dec_proj,decdproj,decmproj,decsproj)
      return
      end



      subroutine hms_to_deg(rah,ram,ras,ra_deg)

      real rah,ram,ras,ra_deg

      ra_deg = 15.*(rah + ram/60. + ras/3600.)

      return
      end



      subroutine dms_to_deg(decd,decm,decs,dec_deg)

      real decd,decm,decs,dec_deg

      if (decd.gt.0.0) then
         dec_deg = decd + decm/60. + decs/3600.
      else
         dec_deg = decd - decm/60. - decs/3600.
      endif

      return
      end



      subroutine deg_to_hms(ra_deg,rah,ram,ras)

      real ra_deg,rah,ram,ras

      rah = aint(ra_deg/15.)
      ram = aint( (ra_deg/15.-rah)*60. )
      ras = ((ra_deg/15. - rah)*60. - ram)*60.

      return
      end




      subroutine deg_to_dms(dec_deg,decd,decm,decs)

      real dec_deg,decd,decm,decs

      sign = 1.0
      if (dec_deg.lt.0.) sign = -1.0
      dummy = sign * dec_deg

      decd = aint(dummy)
      decm = aint((dummy - decd)*60.)
      decs = ((dummy - decd)*60 - decm)*60.
      decd = sign * decd

      return
      end



      subroutine remove_out_of_field_objects(idone)

      include 'maskgen.common.blocks'

      integer i,idone

c		...  check for objects out of field of view: they get thrown out
c		...  regardless of priority, but print warning if one of the stars is
c		... an alignment star.

      if (idone.eq.0) return

      do i=1, icount
         if (projpriority(i).ge.0.0) then
            xedge = endslit
            yedge = globslit/2.0
         else
            xedge= starwidth/2.0
            yedge= starwidth/2.0
         endif
		                  
         if     (xproj(i) .lt. (-0.5*xsize +xedge) 
     +      .or. xproj(i) .gt. (0.5 *xsize -xedge )
     +      .or. yproj(i) .lt. (-0.5*ysize + yedge)
     +      .or. yproj(i) .gt. (0.5* ysize - yedge) ) then

c	      ....      throw out the object if outside the field of view for 
c	      ....      the automated slit selection/algorithm case


           if(auto_priority) then
              ifull = kcount(i)
              if (priority(ifull).lt.0.0 .and. iwarn.ne.1) then
                 write(6,'(41h ****WARNING: alignment star out of field)')
                 write(6,'(15h object number ,i10)')inumobj(ifull)
              endif
              iflag(ifull) =1
              idone =0
           endif
         endif   !   .... done with field of view check for this object

c		@@@@@@@   should examine for wavelength coverage here, and either throw
c			  ... out, or flag with warning, users option... this is an
c			  ... additional test for yproj, but only will apply to slits,
c			  ... not alignment stars...	

      enddo
      return
      end



      subroutine compute_slit_ends(idone)


      include 'maskgen.common.blocks'

      parameter (wfpix = 0.77402)  ! the wfccd pixel size in arcsec

c
c 			.......set up preliminary  slit coordinates.......

c		... NOTE: the variables slitpose,w and slitpose,w proj are
c		... normally zero: they are used for manual adjustment of
c		...slit boundaries: see OPTION 7
c
c		... Set up preliminary set start, stop positions without checking
c		... for endslit or slitlength criteria; however, build in
c		... the box size for the alignment stars

c
c endpix is minimum distance (in pixels) between slits, wfpix is pixel scale
c

        if (idone.eq.0) return

        do i=1,icount
           call get_adjacent_objects(i,ilast,inext)
           if (projpriority(i).ge.0.0) then     ! if i_th object is not an alignment star...

              if ( ilast.ne.0 .and. inext.ne.0 ) then

c this block if adjacent to an alignment box

                 if (projpriority(ilast).lt.0.0) then
                    xstart(i) = xproj(ilast) + starwidth/2. + endpix*wfpix
                 else
                    xstart(i)=(xproj(i)+xproj(ilast))/2.+slitposeproj(i) + endpix*wfpix*0.5
                 endif
                 if (projpriority(inext).lt.0.0) then
                    xend(i) = xproj(inext) - starwidth/2. - endpix*wfpix
                 else
                    xend(i) = (xproj(i)+xproj(inext))/2.+slitposwproj(i) - endpix*wfpix*0.5
                 endif
c
              elseif (ilast .eq. 0 .and. inext .eq. 0) then
	         xstart(i) = -0.5*xsize          ! 1st object slit runs to low edge of field
                 xend(i)   =  0.5*xsize          ! slit runs to high edge of field of view
              elseif (ilast .eq. 0) then
	         xstart(i) = -0.5*xsize       ! 1st object slit runs to low edge of field
                 xend(i) = (xproj(i)+xproj(inext))/2.+slitposwproj(i) - endpix*wfpix*0.5
              elseif (inext .eq. 0) then
                 xend(i) = 0.5*xsize          ! slit runs to high edge of field of view
                 xstart(i)=(xproj(i)+xproj(ilast))/2.+slitposeproj(i) + endpix*wfpix*0.5
              endif
              if (holes) then
                 xstart(i) = xproj(i) - wproj(i)/2.
                 xend(i)   = xproj(i) + wproj(i)/2.
              endif
              call compute_yends(i)
           else                               ! alignment star
              xstart(i) = xproj(i) - starwidth/2.0
              xend(i)   = xproj(i) + starwidth/2.0
              ystart(i) = yproj(i) 
              yend(i)   = yproj(i) 
           endif
c
c 		...get several slit properties...
c               ...note that east and west really mean from  left to right on the plot
c	        ...it is not e to west if rotated, but rather left to right;
c		...the slits always run in the X (left/right) direction...
c
           percent(i) = 1. - (xproj(i)-xstart(i))/(xend(i)-xstart(i))
           lproj(i) = xend(i)-xstart(i)
           lwest(i) = xproj(i)-xstart(i)
           least(i) = xend(i)-xproj(i)


        enddo     !...the slit coordinates are now tentatively set up, but they may
c                  ... violate minimum slit clearance/length criteria, or even overlap.

        return
        end



        subroutine get_adjacent_objects(i,ilast,inext)

c a routine to get the adjacent objects within the stacking area

        include 'maskgen.common.blocks'
        integer i,ilast,inext

        ilast = 0
        inext = 0

c find closest high side object

        xmin = xsize
        do j = i+1,icount
           if (iflag(kcount(j)).ne.1) then
              ydist = abs(yproj(j)-yproj(i)) 
              xdist = abs(xproj(j)-xproj(i)) 
              if (ydist .le. ystack .and. xdist.le.xmin) then
                 inext = j
                 xmin = xdist
              endif
	   endif
	enddo

c now closest low side object

        xmin = xsize
        do j = 1,i-1
           if (iflag(kcount(j)).ne.1) then
              ydist = abs(yproj(j)-yproj(i)) 
              xdist = abs(xproj(j)-xproj(i)) 
              if (ydist .le. ystack .and. xdist.le.xmin) then
                 ilast = j
                 xmin = xdist
              endif
           endif
        enddo

c	endif
	return
	end




      subroutine compute_yends(i)
c
c compute ystart and yend given the cass ring rotation PA, and PA of i-th slit
c
      integer i

      include 'maskgen.common.blocks'
      parameter (raddeg= 57.29577951)

      if (tilted_slits) then
         alpha = (thinrot-180.-paproj(i) )/raddeg
         ystart(i) = yproj(i) - (xproj(i) - xstart(i)) * tan(alpha) 
         yend(i)   = yproj(i) + (xend(i) - xproj(i)) * tan(alpha) 
      else
         ystart(i) = yproj(i)
         yend(i)   = yproj(i)
      endif

      return
      end




      subroutine delete_alignment_clashes(idone)

      include 'maskgen.common.blocks'

      parameter (wfpix = 0.77402)  ! the wfccd pixel size in arcsec

      integer idone,i,j
c
c   if object is an alignment star delete all adjacent objects which clash
c

      if (idone.eq.0) return

      do i=1,icount
         if (projpriority(i).lt.0.0 .and. iflag(kcount(i)).ne.1) then                 
c
c remove objects with more -ve x which clash
c
            dmin = endslit + endpix*wfpix
            do j = i-1,1,-1
               if (iflag(kcount(j)).ne.1) then
                  if ( (xstart(i)-xproj(j)).lt.dmin .and. projpriority(j).ge.0.0 
     %                .and. abs(yproj(i)-yproj(j)).lt.ystack ) then
                     iflag(kcount(j)) = 1
                     idone = 0
                  endif
               endif
            enddo

c
c remove objects with more +ve x which clash
c
            do j = i+1,icount
               if (iflag(kcount(j)).ne.1) then
                  if ( (xproj(j)-xend(i)).lt.dmin .and. projpriority(j).ge.0.0
     %                .and. abs(yproj(i)-yproj(j)).lt.ystack ) then
                     iflag(kcount(j)) = 1
                     idone = 0
                  endif
               endif
            enddo
         endif
      enddo

      return
      end



        subroutine delete_slit_clashes(idone)

        include 'maskgen.common.blocks'

        real dxslit
        integer idone,ifull,i
 
c
c       check the endslit length from the object to the top end of the slit
c
        if (idone.eq.0) return

        do i = 1,icount-1
c
c  check for conflicts unless this slit and next slit are both alignment stars
c
           if (iflag(kcount(i)).ne.1) then
              call get_adjacent_objects(i,ilast,inext)
              if (inext.gt.0)   then
                 if(projpriority(i).ge.0.0 .or. projpriority(inext).ge.0.0) then
                    dxslit = xend(i) - xproj(i)
                    if(dxslit.le.endslit) then         ! not enough dxslit--drop lower priority
                       if(projpriority (i).gt.projpriority(inext)) then
                          ifull = kcount(i)
                       else
                          ifull = kcount(i+1)
                       endif    
                       iflag(ifull) =1   ! dropping an object
                       idone = 0
                    endif
                 endif
              endif
           endif
        enddo  

        return
        end




        subroutine delete_hole_clashes(idone)

        include 'maskgen.common.blocks'

        parameter (wfpix = 0.77402)  ! the wfccd pixel size in arcsec

        real dx
        integer idone,ifull,i
 
        if (idone.eq.0) return

        do i = 1,icount-1
           if (iflag(kcount(i)).ne.1) then
              call get_adjacent_objects(i,ilast,inext)
              if (inext.gt.0)   then
                 if (projpriority(i).ge.0.0 .or. projpriority(i+1).ge.0.0) then
                    dx = xstart(inext) - xend(i)
                    if (dx.lt.endpix*wfpix) then
                       if(projpriority(i).gt.projpriority(inext)) then
                          ifull = kcount(i)
                       else
                          ifull = kcount(inext)
                       endif    
                       iflag(ifull) =1   ! dropping an object
                       idone = 0
                    endif
                 endif
              endif
           endif
        enddo

        return
        end




        subroutine check_slit_lengths(idone)
c
c check that total slit length is at least equal to SLITMIN
c

        include 'maskgen.common.blocks'

        integer idone,i,ifull
        real p0,plo,phi

        if (idone.eq.0) return

        do i = 1,icount
           if (iflag(kcount(i)).ne.1) then
              p0 =projpriority(i)

c               ... check whether this, and adjacent objects are alignboxes or not:

              if (p0.gt.0.0) then  ! this object is not an alignbox...
                 sll = xend(i) - xstart(i)
                 if (sll.lt.slitmin) then  ! this slit length too short
                    call get_adjacent_objects(i,ilast,inext)

c                    --which object should be dropped:

                    if(ilast.eq.0) plo= -1.0
                    if(ilast.ne.0) plo= projpriority(ilast)
                    if(inext.eq.0) phi= -1.0
                    if(inext.ne.0) phi= projpriority(inext)

c               ... remove the lowest priority object of plo,p0,phi (ie highest value of p)

                    pmax = max(plo,p0,phi)
                    if(pmax.eq.plo) ifull = kcount(ilast)
                    if(pmax.eq.p0)  ifull = kcount(i)
                    if(pmax.eq.phi) ifull = kcount(inext)
                    iflag(ifull) =1   ! dropping an object
                    idone=0
                 endif
              endif
           endif
        enddo

        return
        end






        subroutine get_input_file_name(infile)

        character*64 infile

        write(6,'(22h Enter input filename:)')
        read(5,'(a)')infile
        return
        end




        subroutine get_title(title)

        character*32 title

        write(6,'(47h Enter a title (<=32 chars.) for the mask plot:)')
        read(5,'(a32)')title

        return
        end




      subroutine get_coords

      include 'maskgen.common.blocks'

      logical error
      character*80 line

      error = .false.

      do while(.true.)

         error = .false.
         write(6,'(a)')'Enter RA (hms) and Dec (dms) of input center on one line'
         read(5,'(a80)')line
         if (line.eq.'             ') return

c read RA string from line

         i1 = 0
         call get_next_number(i1,i2,rah,line,error)
         i1 = i2
         call get_next_number(i1,i2,ram,line,error)
         i1 = i2
         call get_next_number(i1,i2,ras,line,error)

c read declination string from line

         i1 = i2
         call get_next_number(i1,i2,decd,line,error)
         i1 = i2
         call get_next_number(i1,i2,decm,line,error)
         i1 = i2
         call get_next_number(i1,i2,decs,line,error)

         if (.not.error) then
            write(6,'(a)')'Enter Epoch'
            read(5,*) epoch
            coords = .true.
            return
         endif

      enddo

      return
      end





      subroutine get_next_number(i1,i2,anum,line,error)
c
c find next string delimted by a space or colon and test that it is a number
c
      character*80 line
      logical error,is_a_number

      if (error) return

      ilast = i1
      do ic=i1+1,80
         if (line(ic:ic).eq.' ' .or. line(ic:ic).eq.':') then
            if (ic.gt.ilast+1) then
               if ( is_a_number(line,ilast+1,ic-1) ) then
                  read(line(ilast+1:ic-1),*)anum
                  i2 = ic
               else
                  error = .true.
               endif
               return
            else
               ilast = ic   !if consecutive delimiters change ilast
            endif
         endif
      enddo
      write(6,'(a)')'Error: missing number'
      error = .true.
      return
      end

            

      logical function is_a_number(line,i1,i2) 

      character*80 line
      logical error,not_a_digit

      is_a_number = .true.
      error = .false.

c check for sign in first position
      is = i1
      if (line(i1:i1).eq.'+' .or. line(i1:i1).eq.'-') then
          is = i1+1
          if (is.eq.i2) error = .true.
      endif

c check digits and count decimal places

      ndp = 0
      do j = is,i2
          if (not_a_digit(line(j:j)) .and. line(j:j).ne.'.') error=.true.
          if (line(j:j).eq.'.') ndp = ndp + 1
      enddo

      if (ndp.gt.1) error = .true.

      if (error) then
         write(6,'(a)')'Error: not a numeric string'
         is_a_number = .false.
      endif
      return
      end




      logical function not_a_digit(c)

      character*1 c

      not_a_digit = .true.

      if (c.eq.'1' .or. c.eq.'2' .or. c.eq.'3' .or.
     %    c.eq.'4' .or. c.eq.'5' .or. c.eq.'6' .or.
     %    c.eq.'7' .or. c.eq.'8' .or. c.eq.'9' .or. c.eq.'0') then
          not_a_digit = .false.
      endif
      return
      end


       


      subroutine set_alignment_warning

      include 'maskgen.common.blocks'

      write(6,'(55h Enter 1 to skip warning if alignment star out of field)')
      write(6,'(24h Enter 0 to give warning)')
      read(5,*)iwarn

      return
      end





        subroutine change_defaults

        include 'maskgen.common.blocks'

        parameter (wfpix = 0.77402)  ! the wfccd pixel size in arcsec

        character*1 ans
        integer ians


        write(6,'(51h Enter 99 to modify default parameters, otherwize 0)')
        read (5,*)idef
        if (idef.ne.99) return

        write(6,'(27h Is this mask for COSMIC ? )')
        read (5,*)ans
        if (ans.eq.'y' .or. ans.eq.'Y') then
           cosmic = .true.
           scale = cosmic_scale
        else
           cosmic=.false.
        endif


        if (cosmic.ne..true.) then
           write(6,'(42h Enter 1 to change alignment star box size)')
           write(6,'(45h Enter 0 to accept default size [8.0 arc sec])')
           read(5,*) ians
           if (ians.eq.1) then
              write(6,'(28h Enter alignment star width )')
	      read(5,*)starwidth
           endif
        endif

 
        write(6,'(66h Enter 1 to change normal polarity convention, 0 to accept default)')
	write(6,'(40h default is E is positive, N is positive)')
        read(5,*)isense
        if (isense .eq. 1) then
           write(6,'(42h x axis is normally +ve in the E direction)')
           write(6,'(48h Enter 1 to change sign of x axis in input file:)')
           read(5,*)ixsign
           write(6,'(42h y axis is normally +ve in the N direction)')
           write(6,'(48h Enter 1 to change sign of y axis in input file:)')
           read(5,*)iysign
        endif         


        if (cosmic.ne..true.) then
           write(6,'(46h Enter 1 if input scale in pixels (not arcsec))')
	   write(6,'(48h AND input *NOT* from WFCCD direct image; else 0)')
           read(5,*)iscale
           if (iscale .eq. 1) then

              write(6,'(37h Enter pixel scale (in arcsec/pixel):)')
              read(5,*)platescale

              write(6,'(52h Enter 1 if input in pixels from WFCCD image; else 0)')
              read(5,*)iwfccd
              if (iwfccd.ne.0) 
     %           stop  ' sorry, wfccd input not implemented yet'
           endif
        endif



        write(6,'(52h Enter number of lines to skip at top of input file:)')
        read(5,*)iskip


c		... bypass or set automatic offsets

        write(6,'(51h Enter 0 to bypass automatic offset changes, else 1)')
	read(5,*)ioffset     
        if (ioffset.gt.0) call change_offsets

        write(6,'(40h Enter 0 to set priorities automatically)')
	write(6,'(45h (alignment star squares are always selected))')
	read(5,*)kpri
        if (kpri.eq.0)then
           auto_priority = .true.
        else
           auto_priority = .false.
        endif
c
        if (cosmic.ne..true.) then
           write(6,'(52h Enter 1 to modify default field size limits, else 0)')
           write(6,'(49h default values are 1500 arcs along slits [x] and)')
           write(6,'(41h 1080 arcs along dispersion direction [y])') 
	   read(5,*)imod
	   if (imod.eq.1) then
              write(6,'(39h Enter x,y field size limits in  arcsec)')
	      read(5,*)xsize,ysize
	   endif
        else
           xsize = 710.0
           ysize = 450.0
        endif

        write(6,'(a)')'Enter 0 for slits (default), 1 for holes'
        read(5,*)ihole

c set up hole parameters

        if (ihole.eq.1) then
           holes = .true.
           do while (globslit.lt.1.55 .or. globslit.gt.2.08)
              write(6,'(a)')'Enter hole diamters (1.55 to 2.08 arc sec permitted)'
              read(5,*)globslit
           enddo
           endslit = globslit/2.0
           slitmin = globslit
           write(6,'(a)')'Enter minimum hole edge distance (>= 4.0 pixels recommended)'
           read(5,*)endpix
        else

c set up slit parameters

           write(6,'(a)')'Enter 1 to modify slit separation, slit end-to-object and'
           write(6,'(a)')'minimum slit length'
           write(6,'(a)')'Enter 0 to accept defaults (4.0 pix, 2.0 arcsec, 10.0 arcsec)'
	   read(5,*)imod
           if(imod.eq.1) then
	      write(6,'(62h Enter minimum slit separation (pixels), minimum end-to-object)')
              write(6,'(52h distance (arcsec), and minimum slit length (arcsec))')
              read(5,*)endpix,endslit,slitmin
           endif

           write(6,'(a)')'Are tilted slits required? (requires individual P.A.s in 6th column of input)'
           read(5,'(a1)')ans
           if (ans.eq.'y' .or. ans.eq.'Y') tilted_slits = .true.

        endif

           write(6,'(a)')'Do you wish to stack spectra?'
           read(5,'(a1)')ans
           if (ans.eq.'y' .or. ans.eq.'Y') then
              if (tilted_slits) then
                 write(6,'(a)')'Cannot stack tilted slits'
              else
                 stacked = .true.
                 write(6,'(a)')'Enter minimum y distance between slits in pixels'
                 read(5,*)ystack
                 ystack = ystack * wfpix
              endif
           else
              ystack = ysize
           endif



        return
        end


        subroutine set_global_slitwidth

        include 'maskgen.common.blocks'

        if (holes) return

	write(6,'(37h Enter the global slitwidth in arcsec)')
c       write(6,'(61h (...you can edit individual ones later in the text file only)'
	read(5,*)globslit

        return
        end




        subroutine make_postscript_file(dxpaper,dypaper,numbers)

        include 'maskgen.common.blocks'

        logical numbers
        real tenlarge

c
c plot the current working set of vectors--ioption is initialized to zero
c

        call initialize_photo_plot(dxpaper,dypaper,tenlarge,x2,y2)
        call make_plot_title(tenlarge,x2,y2)
	call plot_slits
	call plot_fiducial_marks(scale)
        if (numbers) then
           call plot_numbers
           numbers = .false.
        endif
        call pgend

        return
        end



        subroutine rename_photomask_pgplot_file

        character*67 photomask,line*80

        write(6,'(a)')'Enter name of photomask postscript file'
	read(5,'(a67)')photomask
        line(:13) = 'mv pgplot.ps '
        line(14:) = photomask
        call system(line)
	return
	end




        subroutine make_xwindow_plot(xsize,ysize,scale)

        real xsize,ysize,scale,enlarge

        call initialize_xwindow_plot(xsize,ysize,scale,enlarge,x1,y1,x2,y2)
        call make_plot_title(enlarge,x2,y2)
	call plot_slits
	call plot_fiducial_marks(scale)
        call plot_numbers
        call pgend

        return
        end




        subroutine initialize_xwindow_plot(xsize,ysize,scale,enlarge,x1,y1,x2,y2)

        real xsize,ysize,scale,tenlarge,x1,y1,x2,y2,facx,facy

        call pgbegin(0,'/xw',1,1)  ! open the xwindow

c       facx = 8.115 * 25.4 * scale / xsize
c       facy = 4.869 * 25.4 * scale / ysize
        facx = 8.50 * 25.4 * scale / xsize
        facy = 6.50 * 25.4 * scale / ysize

        tenlarge = facx
        if (facy.lt.facx) tenlarge = facy

        x1=0.25
        x2=xsize/scale*tenlarge/25.4+0.25
        y1=0.25
        y2=ysize/scale*tenlarge/25.4+0.25

        call pgvsize(x1,x2,y1,y2)
        x1=xsize/(-2.)
        x2=xsize/(2.)
        y1=ysize/(-2.)
        y2=ysize/(2.)
        call pgwindow(x2,x1,y1,y2)
        call pgslw(1)
        call pgbox('bc',0.,0.,'bc',0.,0.)
        return
        end










        subroutine plot_numbers

        include 'maskgen.common.blocks'

        character*5 text

        call pgsch(0.8)
        do i=1,icount
           text='     '
           write(text,'(i5)') iobjproj(i)
           x1=xproj(i)
           y1=yproj(i)+10.
           if (iobjproj(i) .ge. 1000) then
              call pgptext(x1,y1,0.,0.5,text(1:5))
           else if (iobjproj(i) .ge. 100) then
              call pgptext(x1,y1,0.,0.5,text(2:5))
           else if (iobjproj(i) .ge. 10) then
              call pgptext(x1,y1,0.,0.5,text(3:5))
           else if(iobjproj(i).ge.0) then
              call pgptext(x1,y1,0.,0.5,text(4:5))
           else
              call pgptext(x1,y1,0.,0.5,text)
           endif
           y1=yproj(i)
           call pgpoint(1,x1,y1,2)
        enddo
        call pgsch(1.0)
        return
        end





        subroutine initialize_photo_plot(dxpaper,dypaper,tenlarge,x2,y2)

        include 'maskgen.common.blocks'
        real tenlarge

c        ..../ps to create pgplot.ps 
        call pgbegin(0,'/ps',1,1)

        x1 = 0.0
        y1 = 0.02

        xmax = dxpaper - 0.10
        ymax = dypaper - 0.55

        facx = (xmax-x1) * 25.4 * scale / xsize
        facy = (ymax-y1) * 25.4 * scale / ysize
        tenlarge = facy
        if (facx.lt.facy) tenlarge = facx

        x2=xsize/scale*tenlarge/25.4 + x1
        y2=ysize/scale*tenlarge/25.4 + y1
c
        aspect = dypaper/dxpaper
        call pgpaper(dxpaper,aspect) 
        call pgvsize(x1,x2,y1,y2)
        x1=xsize/(-2.)
        x2=xsize/(2.)
        y1=ysize/(-2.)
        y2=ysize/(2.)
        call pgwindow(x2,x1,y1,y2)
        call pgslw(1)
        call pgbox('bc',0.,0.,'bc',0.,0.)

        return
        end








c       subroutine initialize_photo_plot(dxpaper,dypaper,tenlarge,x2,y2)
c
c       include 'maskgen.common.blocks'
c       real tenlarge
c
c        ..../ps to create pgplot.ps 
c       call pgbegin(0,'/ps',1,1)
c
c       dxpaper = 16.5
c       dypaper = 10.5
c       x1 = 0.0
c       y1 = 0.02
c       x2=xsize/scale*tenlarge/25.4 + x1
c       y2=ysize/scale*tenlarge/25.4 + y1
c
c       aspect = dypaper/dxpaper
c       call pgpaper(dxpaper,aspect) 
c       call pgvsize(x1,x2,y1,y2)
c       x1=xsize/(-2.)
c       x2=xsize/(2.)
c       y1=ysize/(-2.)
c       y2=ysize/(2.)
c       call pgwindow(x2,x1,y1,y2)
c       call pgslw(1)
c       call pgbox('bc',0.,0.,'bc',0.,0.)
c
c       return
c       end




        subroutine make_plot_title(tenlarge,x2,y2)

        include 'maskgen.common.blocks'

        real x2,y2
        character*32 string2
        character*40 string3
        character*60 string4

c 
c put title and useful text in top left corner
c

        call pgsch(0.6)
c       call pgtext(x2-20.,y2+42.0,title)   ! -12.0 to +48.0
        call pgtext(x2-20.,y2+47.0,title)   ! -12.0 to +48.0

        write(string2,'(14hMagnification=,f5.2)') tenlarge
c       call pgtext(x2-20.,y2+30.0,string2)   ! -24.0 to +36.0
        call pgtext(x2-20.,y2+33.0,string2)   ! -24.0 to +36.0

 	if(thetad.gt.359.99) thetad = thetad- 360.0
        write(string3,155) thinrot, thetad
 155    format('WFCCD PA = ',f7.2,' Input PA= ',f7.2)
c       call pgtext(x2-20.,y2+18.0,string3)    ! -36.0 to +24.0
        call pgtext(x2-20.,y2+19.0,string3)    ! -36.0 to +24.0

        if (coords) then
           irah = nint(rahproj)
           iram = nint(ramproj)
           idecd = nint(decdproj)
           idecm = nint(decmproj)

           dumdecs = decsproj
           if ( abs(60.0-decsproj).lt.0.05 ) then
              dumdecs = 0.0
              idecm = idecm + 1
           endif
           dumras = rasproj
           if ( abs(60.0-rasproj).lt.0.005 ) then
              dumras = 0.0
              iram = iram + 1
           endif

           write(string4,150)irah,iram,dumras,idecd,idecm,dumdecs,epoch
        else
           write(string4,152) x0unproj,y0unproj
        endif
 150    format('RA ',i2,x,i2,x,f5.2,5x,'DEC 'i3,x,i2,x,f4.1,x,
     %         'EPOCH ',f6.1) 
 152    format('x0=',f7.2,'" y0=',f7.2,'" relative to input (x=0,y=0)')
c       call pgtext(x2-20.,y2+6.0,string4)  ! -48.0 to +12.0
        call pgtext(x2-20.,y2+5.0,string4)  ! -48.0 to +12.0
        call pgsch(1.0)

        return
        end






	subroutine plot_slits

        include 'maskgen.common.blocks'


        real xpts(5),ypts(5),w

        do i=1,icount

c for tilted slits use the projected y-width

           if (holes .and. projpriority(i).ge.0.0) then
              call plot_hole(xproj(i),yproj(i),wproj(i))
           else
              w = wproj(i)
              if (ystart(i).ne.yend(i)) then
                 w=w*sqrt(((ystart(i)-yend(i))/(xstart(i)-xend(i)))**2+1.0)
              endif
           
              xpts(1)=xstart(i)
              ypts(1)=ystart(i) - w/2.

              xpts(2)=xstart(i)
              ypts(2)=ystart(i) + w/2.

              xpts(3)=xend(i)
              ypts(3)=yend(i) + w/2.

              xpts(4)=xend(i)
              ypts(4)=yend(i) - w/2.

              xpts(5)=xstart(i)
              ypts(5)=ystart(i) - w/2.

c		... number of line strokes
              call pgslw(1)
              call pgsfs(1)
              call pgpoly(5,xpts,ypts)
              call pgslw(1)
           endif
        enddo
        return
        end




        subroutine plot_hole(xc,yc,w)

c plot a 24-sided polygon, radius w/2 to represent hole

        parameter (raddeg= 57.29577951)
        real x(25),y(25),xc,yc,w,r

        r = w/2.
        do i = 1,25
           alpha = real(i-1)*360./24.
           x(i) = xc + r*cos(alpha/raddeg)
           y(i) = yc + r*sin(alpha/raddeg)
        enddo
c       
        call pgslw(1)
        call pgsfs(1)
        call pgpoly(25,x,y)
        call pgslw(1)
        return
        end




	subroutine plot_fiducial_marks(scale)

        real xgrid(4),ygrid(4),scale

c
c       Put 4 fidicucial marks on: 2 pairs spaced 300 mm on the output in x, and
c       10 mm above southern boundary (which is negative y)
c       the separation in y should average 60 mm times scale
c       the separation in x should be 80mm times scale.
cc      225 mm above the first cross in Y:  photo mask distances should be
cc      80mm and 60 mm separation



        call pgsch(2.0)  ! make the cross a little bigger
        xgrid(1)= - 40.0*scale
        ygrid(1)= -330.0  ! in arcsec

        xgrid(2)= xgrid(1)
        ygrid(2)=  ygrid(1) + 60.0*scale


        xgrid(3)= - xgrid(1)
        ygrid(3)= ygrid(1)

        xgrid(4) = xgrid(3)
        ygrid(4) = ygrid(2)
        call pgpoint(4,xgrid,ygrid,2)
        call pgsch(1.0) ! return it to its old size

        return
        end


        subroutine get_commands(ioption)

        integer ioption

        ioption = 0
        write(6,'(59h Enter 1 to change the rotation angle and re-select objects)')
        write(6,'(29h      2 to exclude an object,)')
        write(6,'(29h      3 to include an object,)')
        write(6,'(38h      4 to plot again on the terminal,)')
        write(6,'(45h      5 type current list to term/outputfile,)')
        write(6,'(39h      6 to write a photomask.ps file   )')
        write(6,'(32h      7 to edit lengths of slits)')
        write(6,'(30h      8 to include all objects)')
        write(6,'(30h      9 to exclude all objects)')
        write(6,'(53h      10 to add overall offset in x or y to all slits,
     %        14h and reselect:)')
        write(6,'(44h      11 to exit program, no further output:,//)')
        read(5,*)ioption

        return
        end




        subroutine change_cass_ring_pa

        include 'maskgen.common.blocks'

        parameter (raddeg= 57.29577951)

        do while (.true.)
           write(6,'(48h Enter the CASS RING angle for your observations)')
	   write(6,'(48h acceptable angles are 90.0 < ringangle <= 270.0)')
	   write(6,'(51h 180 gives N/S slits; 270 recommended for E/W slits)')
 180       format('Centroid x0=',f8.2,'" and y0=',f8.2,'" relative to the',
     +     ' input (x=0,y=0) coordinates; Rotation angle:  E=',f8.2,
     +     ' degree CW from positive x-axis','  Set COSMIC base angle',
     +     ' to ',f8.2)
	   read(5,*)thinrot
           if (thinrot.ge.89.9999 .and. thinrot.lt.270.0001) then
	      thetad = thinrot+90.0
	      theta  = thetad/raddeg
              return
           endif
	   write(6,'(32h illegal CASS ring angle entered)')
	enddo

        return
        end




        subroutine print_current_list
c
c prints object information to terminal, and to a file if required
c
        logical hardcopy_requested
        integer iunit

        iunit=6
        call print_list(iunit)
        if (hardcopy_requested()) then
           iunit=31
           call print_list(iunit)
           close(unit=31)
        endif
        return
        end



        subroutine print_list(iunit)
c
c print object information to iunit
c
        include 'maskgen.common.blocks'

        integer iunit

        write(iunit,'(a)') ' Current Object Status'
        write(iunit,*) icount
        write(iunit,'(a)') 'Included objects:'
        write(iunit,'(a)')
     %     'Object     X        Y        Width   length   percent  L(West)  L(East)'
        write(iunit,'(a)')
     %     '         (arcs)   (arcs)    (arcs)   (arcs)             (arcs)   (arcs)'
        write(iunit,'(a)')  '  '

        do j = 1,icount
           write(iunit,'(i6,7f9.2)') iobjproj(j),xproj(j),
     %     yproj(j),wproj(j),
     +     lproj(j),percent(j),lwest(j),least(j)
        enddo
        write(iunit,'(a)')  '  '
        write(iunit,'(a)') 'Excluded objects (unrotated postions): '
        write(iunit,'(a)') 'Object   X(arcs)  Y(arcs)    SW(arcs)   PA'

        write(iunit,'(a)') '  '
        do i=1,inumin
           if (iflag(i) .eq. 1) then
              write(iunit,104) inumobj(i),xobj(i),yobj(i),wslit(i),pa(i)
           endif
        enddo
 104    format(i5,x,4(f9.2,x))
        return
        end
 


        logical function hardcopy_requested()
c
c asks user if hardcopy of object info is required in a local file.
c file is opened at unit=31
c
        character*64 tempfile
        integer ktemp

        hardcopy_requested = .false.
        
        write(6,'(61h enter 1 to also write to an output file (not the final file))')
	read(5,*)ktemp

        if (ktemp.eq.1) then    
           hardcopy_requested = .true.
           write(6,'(a)')' Enter file name for current object info'
           read(5,'(a)')tempfile
           open(unit=31,file=tempfile)
        endif
        return
        end



      subroutine select_all_objects(iflag,inumin)

      parameter (kmax=4000)  ! maximum number of objects in a mask list

      integer iflag(kmax),inumin

c 			flag=0 include; flag=1 exclude

      do i=1,inumin
         iflag(i)=0
      enddo
      return
      end



        subroutine exclude_all_objects(iflag,inumin)

        parameter (kmax=4000)  ! maximum number of objects in a mask list

        integer iflag(kmax),inumin

c 			flag=0 include; flag=1 exclude

        do i=1,inumin
           iflag(i)=1
        enddo
        return
        end




        subroutine change_offsets

        include 'maskgen.common.blocks'

        integer ioffset
        real xoff,yoff

     
        write(6,'(a)')'Select offset options:'
        write(6,'(a)')'0: no offset,no centroid'
        write(6,'(a)')'1: offset only'
        write(6,'(a)')'2: only centroid'
        write(6,'(a)')'3: offset and centroid'
        read(5,*)ioffset
	if(ioffset.eq.1.or.ioffset.eq.3) then
		write(6,'(a)')' Enter xoff, yoff in arcsec'
		read(5,*)xoff,yoff
	endif
        return
        end




        subroutine move_slit_boundaries(inumobj,inumin,slitposw,slitpose)

        parameter (kmax=4000)  ! maximum number of objects in a mask list

        real slitposw(kmax),slitpose(kmax),interchange
        integer inumobj(kmax),inumin,iobjw,iobje

        write(6,'(a)')'Enter object number on RIGHT side of slit interface:'
        read(5,*)iobjw
        write(6,'(a)')'Enter object number on LEFT side of slit interface:'
        read(5,*)iobje
        write(6,'(a)')'Enter amount to change interface (+ to move right):'
        read(5,*)interchange
        do i=1,inumin
           if (inumobj(i) .eq. iobjw) then
              slitposw(i)=interchange
           endif
        enddo
        do i=1,inumin
           if (inumobj(i) .eq. iobje) then
              slitpose(i)=interchange
           endif
        enddo
        return
        end



        subroutine print_excluded_objects
c
c    print out list of objects not on mask
c
        include 'maskgen.common.blocks'

        character*64 newfile

        write(6,'(a)')' write input list of all objects not on current mask'
        write(6,'(a)')' but including alignment stars: [1=y, 0=n]'
        read(5,*)ilist
        if(ilist.eq.1) then
           write(6,'(a)')' Enter filename for new list'
           read(5,'(a)')newfile
           open(unit=40,file=newfile,status='unknown')
           jnew=0
           write(40,'(59hExcluded objects and alignment stars (unrotated positions) )')
           write(40,'(a32)')title
           write(40,'(a)')'   ID   X(arsc)  Y(arsc)   Priority   Mag.    P.A.'
           do i=1,inumin
              if(iflag(i).eq.1.or.priority(i).lt.0.0) then
                 jnew = jnew + 1
                 write(40,'(i6,5f9.2)') inumobj(i),xobj(i),
     %                    yobj(i),priority(i),xmag(i),pa(i)
              endif
           enddo
           write(6,'(a,i6)')' # of objects in new list: ', jnew
        endif
        return
        end



        subroutine prompt_for_numbers(numbers)

        logical numbers
        character*1 ans
 
c
c old version...
c       write(6,'(a)')'Enter 1 to include object numbers in plot (for reference purposes)'
c       write(6,'(a)')'Enter 2 for no numbers (for the actual mask):'
c       read(5,*)inumbers
c

        numbers = .false.
        write(6,'(a)')'Do you want numbers plotted on the postscript file?'
        read(5,*)ans
        if (ans.eq.'y' .or. ans.eq.'Y') numbers = .true.
        return
        end




        subroutine produce_gmcode_file

        include 'maskgen.common.blocks'

        real xs(kmax),ys(kmax),xe(kmax),ye(kmax),sw(kmax)
c
c Put slit dimensions in inches, transform to GMC x-direction and 
c order with increasing x(GMC)
c
c       scalein = 3.60671e-03   ! inches per arcsec
        scalein = 3.606689e-03   ! inches per arcsec
        if (cosmic) scalein = 3.33362e-03    ! scale for cosmic

        do i = 1,icount
           j = icount - i + 1
           xs(i) = -xstart(j) * scalein
           ys(i) =  ystart(j) * scalein
           xe(i) = -xend(j) * scalein
           ye(i) =  yend(j) * scalein
           sw(i) =  wproj(j) * scalein
c
c
c check for slit falling off edge of mask
c
           dxdy = (xe(i)-xs(i))/(ye(i)-ys(i))
 
           if (ye(i).gt.0.5*ysize*scalein) then
              
              dye   =  ye(i) - 0.5 * ysize * scalein
              xe(i) = xe(i) - dye * dxdy
              ye(i) = 0.5*ysize * scalein
           endif


           if (ys(i).lt.-0.5*ysize*scalein) then
              dys   = ys(i) + 0.5*ysize * scalein 
              xs(i) = xs(i) + dys * dxdy
              ys(i) = - 0.5*ysize * scalein 
           endif
           
        enddo

	call gmcode_2p(xs,ys,xe,ye,sw,icount,title)

	return
	end






	subroutine create_slit_mask
c
c compute slit positions, apply shifts and check for clashes and bad slits
c and iterate until no more objects are thrown out
c

        include 'maskgen.common.blocks'
        integer idone,iterate

c
c    when idone= 1 we are done, if still throwing out objects idone =0 
c
        idone = 0
        iterate = 0
        do while(idone.eq.0)
           iterate = iterate +1
           idone = 1
           call make_new_xylist
           write(6,110)iterate, icount, x0unproj, y0unproj     ! is this really needed???
           call sort_slits
           call shift_xy                                       ! apply x,y offsets
           call remove_out_of_field_objects(idone)
           call compute_slit_ends(idone)
           call delete_alignment_clashes(idone)

           if (auto_priority) then
              if (holes) then
                 call delete_hole_clashes(idone)
              else
                 call delete_slit_clashes(idone)
                 call check_slit_lengths(idone)
              endif
           endif
        enddo

        write(6,'(37h Number of objects in current mask = ,i5)')icount
 110    format('iteration #, current # of obj, Centroid x,y',2i5, 2f7.2)
	return
	end





        subroutine produce_output_files

        include 'maskgen.common.blocks'

        real dxpaper,dypaper
        logical numbers

c
c identification postscript chart
c
        numbers = .true.
        dxpaper = 10.5
        dypaper =  8.0

        call make_postscript_file(dxpaper,dypaper,numbers)
        call rename_pgplot_file
        call produce_gmcode_file
        call print_excluded_objects
	return
	end


        subroutine rename_pgplot_file

        character*67 idfile,line*80

        write(6,'(a)')'Enter name of ID postscript file'
	read(5,'(a67)')idfile
        line(:13) = 'mv pgplot.ps '
        line(14:) = idfile
        call system(line)
	return
	end



C************************************************************************

	SUBROUTINE gmcode_2p(xs,ys,xe,ye,sw,NS,title)

c	written by Barbara L. Weymann and Ray Weymann

c	this file in: /home/mozart/rjw/rjwadm/ociw/fpcam/cur_softw/slits/gmcode_2p.f
c	revision: Feb.  13th 1999
c	this code is being modified to make two single passes: one with the
c	0.020 or 0.025  end mill and the second pass with the engraving tool. 
C	LAST REVISION Feb. 13, 1999
c       modified by S. C. Trager, 3/3/99, to include tilted slits
c       modified by A.McWilliam 18,March,1999 for subroutine version and tilted slits



C  This program creates an ASCII text file in the code
C  used for the computer that runs the machine that cuts
C  rectangular slits (including
c  alignment star squares, for the masks.
c	          ... Zero point located
c			... at the TOP SURFACE of slitmask blank
c			... at the center of the mask.

c	INPUTS FROM THE TERMINAL:
c       1) The name of the output file:
c	
c       2) The user may change end-mill diameter and default distances.
c		
c	Defaults are: metal blank thickness 0.016, engraving tool tip 0.003
c		angle is  2:1

c	All these dimensions are in INCHES. If sl <= 0.001 a round aperture
c	... of diameter sw is made, at position xs, ys. If sl > 0.001 a 
c	... slit is made whose STARTING position is xs, ys and whose ENDING
c       ... position is xe, ye.  This allows for tilted slits.
c	... origin 0.0, 0.0 is at center of mask
c	...>>>>>> NOTE NOTE NOTE: these x co-coordinates increase to the RIGHT
c	...and these y coordinates increase AWAY from the operater as the mask
c	...is mounted 		

c	INTEGER LINE,NPOINTS,NSWEEPS,NZIGZAGS,I,IFILE  
 
        parameter(kmax=4000)
	REAL*4  XS(kmax),YS(kmax),XE(kmax),YE(kmax),SW(kmax)

	CHARACTER ANS*1
	CHARACTER INFILE*32, OUTFILE1*12,OUTFILE2*12, OUTFL*8
	CHARACTER title*32
c
c default milling parameters
c
	ZMETAL = 0.016
	TIP   =  0.003
	ENDMILLD = 0.020
	UNCUT  = 0.004
c       ZSAFE  = 0.100
        ZSAFE  = 0.060
	ZSTART = 0.01    !may want to change zstart to 0.005

	write(6,'(a)')'Do you wish to check/alter the default milling parameters?'
	read(5,'(a1)')ans

	if (ans.eq.'y' .or. ans.eq.'Y') then
c					....metal mask thickness....
	   WRITE(6,'(a)')'ACCEPT DEFAULT MASK THICKNESS OF 0.016 INCHES? [Y/N]' 
	   READ(5,'(A1)') ANS
	   IF (ANS.EQ.'N'.OR.ANS.EQ.'n') THEN
	      WRITE(6,'(a)')'ENTER THICKNESS IN INCHES'
	      READ(5,*) ZMETAL
	   ENDIF

	   WRITE(6,'(A)')'ACCEPT DEFAULT TOOL TIP OF 0.003 INCHES? [Y/N]' 
	   READ *, ANS
	   IF (ANS.EQ.'N'.OR.ANS.EQ.'n') THEN
	      WRITE(6,'(A)')'ENTER TIP IN INCHES'
	      READ(5,*) TIP
	   ENDIF
	
	   WRITE(6,'(A)')'ACCEPT DEFAULT uncut endmill depth of 0.004"? [Y/N]' 
	   READ *, ANS
	   IF (ANS.EQ.'N'.OR.ANS.EQ.'n') THEN
	      WRITE(6,'(A)')'uncut endmill depth in inches'
	      READ(5,*) UNCUT
	   ENDIF

c          this is the depth below the surface
c          ... to which the endmill cuts.
c          ... this is a POSITIVE number

	endif

	ZENDMILL = ZMETAL - UNCUT 

C       ...Open output file

	WRITE(6,'(A)')' ENTER OUTPUT FILENAME FOR GM-CODE'
	WRITE(6,'(A)')' USE 8 OR FEWER ALPHA-NUMERIC CHARACTERS'
	READ(5,'(A)') OUTFL
	do j=1,8
	   if (outfl(j:j).eq.' ') goto 1
	enddo
 1	len=min(j-1,8)
	OUTFILE1 = OUTFL(1:len) // '.ps1' ! endmill file
	OPEN (UNIT=30,FILE=OUTFILE1,STATUS='UNKNOWN')
	OUTFILE2 = OUTFL(1:len) // '.ps2' ! engraving tool file
	OPEN (UNIT=20,FILE=OUTFILE2,STATUS='UNKNOWN')

	WRITE(20,'(a,a)')   '(', title
	WRITE(20,'(a)') '( 2nd pass: ENGRAVING TOOL'
	WRITE (20, '(A,A)') '(GMFILE: ',OUTFILE2
        WRITE (20, '(A,A)') '(INPUT FILE: ', INFILE
	WRITE(30,'(a,a)')   '(', title
	WRITE(30,'(a)') '( 1st pass: ENDMILL'
	WRITE (30, '(A,A)') '(GMFILE: ',OUTFILE1
        WRITE (30, '(A,A)') '(INPUT FILE: ', INFILE



c	...initialize ...for engraving tool
	WRITE(20,'(a)') 'G4 F 15'    ! pause for 15 seconds before 1st commands...
	WRITE(20,'(a)') 'G98 G70'
	WRITE(20,'(a)') 'G90'
	WRITE(20,'(a)') 'G00'
	WRITE(20,'(a)') 'M03 M07'
	WRITE(20,'(a)') 'F1 E0.1'  ! note the Z motion is 10% of the X-Y motion
	WRITE(20,'(a,f5.1)') ' Z ', ZSAFE
	WRITE(20,'(a)') 'G39 X -2.3 Y - 1.7 I 4.6 J 3.4' ! this displays the full slit
c							   ....pattern

c	...initialize ...for endmill
	WRITE(30,'(a)') 'G4 F 15'    ! pause for 15 seconds before 1st commands...
	WRITE(30,'(a)') 'G98 G70'
	WRITE(30,'(a)') 'G90'
	WRITE(30,'(a)') 'G00'
	WRITE(30,'(a)') 'M03 M07'
	WRITE(30,'(a)') 'F1 E0.1'  ! note the Z motion is 10% of the X-Y motion
	WRITE(30,'(a,f5.1)') ' Z ', ZSAFE
	WRITE(30,'(a)') 'G39 X -2.3 Y - 1.7 I 4.6 J 3.4' ! this displays the full slit
c							   ....pattern

      

c		...Read in array of slits parameters
c		...these should be ordered in
c		... order of increasing X
C 		...XS,YS are starting points at the left of the slit,
C 		...XE,YE are ending points at the right of the slit,
c               ... SW is the slit width, perpendicular to the direction of the slit (AMcW);
c		... x increases to the right. Note that X(servo) increase to RIGHT, whereas
c		... the photographic mask coords increase to LEFT. See the code MAKEMET.F
c
c removed input commands: AMcW 9/march/99
c
	IF(NS.GT.kmax) STOP ' TOO MANY SLITS'
	DO MS =1,2		! first pass is for endmill,
	   DO KS =1, NS
              CALL SLIT(XS(KS),YS(KS),XE(KS),YE(KS),SW(KS),ZSAFE,
     %                  ZSTART,ZMETAL,TIP,MS,ZENDMILL)
	   ENDDO
	ENDDO
	
c	
	WRITE(20,'(a)') ' Z 0.80'
	WRITE(20,'(a)') ' X 2.5  Y -1.6'

c       ... these last two commands move the tool out of the way,
c       ... but in practice against the limits..after the 2nd pass



	
c       ... this is the park position after the 1st pass	
	WRITE(30,'(a)') ' Z 0.80'
	WRITE(30,'(a)') ' X 0.0  Y -1.6' ! as per Oscar request




	WRITE(20,'(a)') ' M30'
	WRITE(30,'(a)') ' M30'

	RETURN
	END
C*******************************************
c       modified 3/3/99 by SCT to allow for tilted slits
c 	... this is the single new routine for
c	... the metal mask GMCODE... it incorporates 

	SUBROUTINE SLIT(XS,YS,XE,YE,SW,ZSAFE,ZSTART,ZMETAL,TIP,
     %              MS, ZENDMILL)


c	... MS=1: This is the 1st (endmill) pass. It writes to file 30
c	... MS=2: This is the 2nd (engraving tool) pass. It writes to file 20

c	Cuts a rectangular slit using the following sequence:
c	1) Sets ABSOLUTE and RAPID mode, and
c	positions the  the tool tip at absolute height ZSAFE above the
c	surface of the blank 
c	---->>>> note these two commands are redundant since they are also
c	the last commands in the subroutine: They are left in to insure that
c	any fiddling outside the routine doesnt cause the machine to 
c	be in the wrong mode or unsafe position. 
c	redundant for safety, since these are also the last commands
c	exiting the routine

C	... If the  SLIT LENGTH , SL is le 0.000 then a hole is made,
c	... but large holes which would require the tool depth to be
c	... under ZUNDERMAX are not allowed.

	character XC*3, YC*3,ZC*3
	character G00*4, G01*4,G90*4,G91*4
	parameter(XC = 'X  ',YC= 'Y  ', ZC = 'Z  ')
	parameter(G00 = 'G00 ', G01 = 'G01 ')
	parameter(G90 = 'G90 ', G91 = 'G91 ')

	ZUNDERMAX = 0.009004   !the maximum tool depth below mask for "narrow" mode
c		... the value of ZUNDERMAX may have to be reduced to 0.008
	ZUNDERNOM = 0.003    !the nominal tool depth below mask for "wide" mode


c switch coordinate order if xe < xs      AMcW 3/99

        IF (XE.LT.XS) THEN
           XDUM = XS
           XS   = XE
           XE   = XDUM
           YDUM = YS
           YS = YE
           YE = YDUM
        ENDIF


	IF(MS.EQ.1) THEN	! ...................the endmill pass starts here .............
	   
	   
	   WRITE(30,'(a4)') G00	!  ensure in RAPID
	   WRITE(30,'(a4)') G90	!  ensure in ABSOLUTE MODE
	   WRITE(30,'(a4,f9.5)') ZC,ZSAFE !  ensure at safe height
	   WRITE(30,'(2(a4,f9.5))') XC,XS,YC,YS ! move to X,Y start position
	   WRITE(30,'(a4,f9.5)') ZC,ZSTART ! move to Z start position rapidly
	   WRITE(30,'(a4)') G01	!  switch to LINEAR (slow) mode
	   WRITE(30,'(a4)') G91	!  switch to INCREMENTAL MODE
		 
	   ZENDINC= -(ZSTART +ZENDMILL)	! the incremental zmotion for endmill
	   WRITE(30,'(a4,f9.5)') ZC,ZENDINC 

 101	   format(2(a4,f9.5))
 102	   format(a4,f9.5)

	   IF(SW.le.0.010) then ! slit is narrow enough to simply make a single cut
c       ... with no small dx, dy offsets to the depth
c       ! cut to X,Y end position
	      WRITE(30,101) XC,(XE-XS),YC,(YE-YS)
	   ELSE IF(sw.gt.0.010) then
c       ... large slit requires moving in a rectangular pattern... 

c  move in Y by projected slit width for tilted slits AMcW 3/99

              SWP = SW * SQRT( ((YE-YS)/(XE-XS))**2 + 1.0 )
              HSWP = 0.5*SWP
c             HSW = 0.5*SW         ! ...half the slit width deleted AMcW 3/99
	      WRITE(30,102) YC,HSWP ! from start to top left
	      WRITE(30,101) XC,(XE-XS),YC,(YE-YS) ! from top left to top right
	      WRITE(30,102) YC,-SWP ! from top right to bottom right
C       from bottom right to bottom left
	      WRITE(30,101) XC,-(XE-XS),YC,-(YE-YS)
	      WRITE(30,102) YC,HSWP ! from bottom left back to start 
	   ENDIF		! end if for narrow vs. wide
	   WRITE(30,'(a4)') G00	!  change to RAPID MODE
	   WRITE(30,'(a4)') G90	!  change to ABSOLUTE MODE
	   WRITE(30,102) ZC,ZSAFE !  move to safe height


	ELSE IF(MS.EQ.2) THEN	!  the 2nd pass with the engraving tool

c	... first calculate the incremental moves and decide between
c	... "narrow" or "wide" case:

c	... calculate the depth below the surface assuming narrow mode:
	   ZUNDER = 2*(SW - TIP) 
	   IF(ZUNDER.LE.ZUNDERMAX) THEN ! use the "narrow mode"
	      KNW =1		! narrow mode
c       ... ZUNDER is a positive #
	   ELSE IF(ZUNDER.GT.ZUNDERMAX) THEN ! use the "wide mode"
	      KNW =2		! wide mode
	      ZUNDER = ZUNDERNOM
	   ENDIF
	   RAD    = 0.5*TIP + 0.25*ZUNDER  
c       ... the radius of the cutting edge when tip is @ ZUNDER
	   DX =     RAD		! shift the start of the cut so that edge of
c       cut is at XS
           XL   = XE-XS
 	   XCUT = XL -2.0*DX	! total length of cut in x direction so that
c       ... total slit length is ~ XL
	   WRITE(20,'(a4)') G00	!  ensure in RAPID
	   WRITE(20,'(a4)') G90	!  ensure in ABSOLUTE MODE
	   WRITE(20,102) ZC,ZSAFE !  ensure at safe height
           WRITE(20,101) XC,XS,YC,YS ! move to X,Y start position
	   
	   ZSTART2 = ZSTART - ZENDMILL ! this is the absolute coordinate
c       ....  starting position for the
c       ... engraving mill pass, which is lower than
c       ... ZSTART by the amount by
c       ... which the endmill has taken out
c       ... the sides of the
c       ... engraving tool should clear
	   ZBELOW = -(ZMETAL + ZUNDER) !... the signed value at the bottom.
	   ZINC   =  ZBELOW - ZSTART2 !... the (negative) incremental Z motion
	   
	   WRITE(20,'(a4,f9.5)') ZC,ZSTART2 ! move to Z start position rapidly

	   WRITE(20,'(a4)') G01	!  switch to LINEAR (slow) mode
	   WRITE(20,'(a4)') G91	!  switch to INCREMENTAL MODE

           DDY = (YE-YS)*DX/XL  !  change in Y position due to DX shift (to keep PA)


	   IF(KNW.EQ.1) THEN	!        narrow mode

	      WRITE(20,101) XC,DX,YC,DDY ! move off start position by DX and DDY
	      WRITE(20,102) ZC,ZINC !  cut  to final zdepth
c             IF (XCUT.NE.0.0 .AND. (YE-YS-2.*DDY).NE.0.0) THEN   ! if not a hole

              IF ( XCUT**2+(YE-YS-2.*DDY)**2 .NE. 0.0 ) THEN   ! if not a hole
	         WRITE(20,101) XC,XCUT,YC,(YE-YS-2.*DDY) !  cut to end of slit
              ENDIF

	   ELSE IF(KNW.EQ.2) THEN !  box mode

c       the total increment in Y necessary to get the desired
c       slit width
              SWP = SW * SQRT( ((YE-YS)/(XE-XS))**2 + 1.0 )
	      DY = SWP - 2*RAD	
c             DY = SW - 2*RAD    ! replaced with above code AMcW 3/99
	      IF(DY.LT.0)
     1          STOP ' negative DY should not occur in subroutine SLIT'
	      DY2 = 0.5*DY	 ! half DY

	      WRITE(20,102) XC,DX   ! move off start position by DX 
	      WRITE(20,102) ZC,ZINC !  cut  to final zdepth
	      WRITE(20,102) YC,DY2 ! cut from center to top left
c       cut from top left to top right
	      WRITE(20,101) XC,XCUT,YC,(YE-YS-2.*DDY)
	      WRITE(20,102) YC,-DY !    cut from top right to bottom right
c       cut from bottom right to bottom left
	      WRITE(20,101) XC,-XCUT,YC,-(YE-YS-2.*DDY)
	      WRITE(20,102) YC,DY2 !    cut from bottom left to start
	   ENDIF		! endif for narrow/wide engraving tool
c	... return to ABSOLUTE and RAPID move
c	... and lift tool to safe Z height

	   WRITE(20,'(a4)') G00	!  change to RAPID MODE
	   WRITE(20,'(a4)') G90	!  change to ABSOLUTE MODE
	   WRITE(20,102) ZC,ZSAFE !  move to safe height
	   
	ENDIF			! ms enddo

	RETURN
	END
C************************************************************************