program cosmicslitmask
c
c program to create a 4x/5x magnification of the negative of a slit mask
c     for use with COSMIC on P200
c the output of this program on a laserprinter will then be photoreduced
c     by 4x/5x in order to create a film (slits clear; everything else
c     dark) that is the mask itself
c
c written by:
c     Mike Pahre
c     Caltech
c     13 February 1995
c
c compile by having the pgplot library in your path; this can be accomplished
c     at Caltech by including the lines in your .cshrc file :
c
c set path = ( $path /usr/local/vlb/pgplot /usr/local/bin/X11 )
c setenv LD_LIBRARY_PATH /usr/lang/SC0.0:/usr/openwin/lib:/usr/local/vlb/pgplot
c
c then the compilation is done by:
c     f77 -O -o cosmicslitmask cosmicslitmask.o -lpgplot -lX11
c
c the program requires an input file with the following specifications:
c   (for a free-format read)
c     col(1)  integer number of object (not necessarily starting at 1)
c     col(2)  real xposition/arcsec; E is positive
c     col(3)  real yposition/arcsec; N is positive
c     col(4)  real width of slit/arcsec
c     col(5)  real magnitude of object
c
c the program will output:
c     (1) the plot of the slitmask; this goes to the terminal if
c         /xterm is chosen; it goes into pgplot.ps if /vps is chosen
c     (2) an output table with the chosen slitmask positions 
c         col(1)  number of object
c         col(2)  xposition/arcsec; E is positive (rot angle=0)
c         col(3)  yposition/arcsec; N is positive (rot angle=0)
c         col(4)  width of slit/arcsec
c         col(5)  length of slit/arcsec
c         col(6)  percentage of the length of slit to find object
c                   0% -> top end of slit; 100% -> S end of slit
c         col(7)  xposition/arcsec relative to centroid (using theta)
c         col(8)  yposition/arcsec relative to centroid (using theta)
c         col(9)  length of slit/arcsec west of object
c         col(10) length of slit/arcsec east of object
c
c note that the output of this program can also be used as the input of
c     a successive iteration of the program itself
c

      integer k
      parameter (k=10000)
      integer i,inumin,inumnow,inumobj(k),iobjproj(k),iflag(k)
      real x0,y0,xobj(k),yobj(k),wslit(k),scale,enlarge
      real xstart(k),xend(k),xproj(k),yproj(k),pi,wproj(k)
      real ystart(k),yend(k),percent(k),lproj(k),xsize,ysize,yoff
      real x1,x2,y1,y2,xpts(5),ypts(5),platescale,interchange,xoff
      real slitposw(k),slitpose(k),slitposwproj(k),slitposeproj(k)
      real least(k),lwest(k),x0unproj,y0unproj,magin(k),magproj(k)
      parameter (scale=11.90)
      parameter (pi=3.1415926535)
      character*64 infile,outfile
      character*32 string1,string2
      character*40 string3
      character*60 string4
      character*4 text

      ioption=0
      inumbers=2
      xoff=0.
      yoff=0.

c
c scale is 11.90 "/mm at prime focus (measured)
c enlarge is enlargement factor for laserprinter plot = 4.0 usually
c

c
c query input filename from user
c
      type *, 'Enter input filename:'
      accept *, infile
      type *, 'Enter a rotation angle (CCW from input x,y; degrees):'
      accept *, theta
      theta=theta/360.*2.*pi
      type *, 'Enter the name for the plot:'
      accept *, string1
      type *, 'Enter 1 to change expected sense of input coordinates:'
      accept *, isense
      if (isense .eq. 1) then
         type *, 'x axis is normally +ve in the E direction;'
         type *, '  Enter 1 to change sign of x axis in input file:'
         accept *, ixsign
         type *, 'y axis is normally +ve in the N direction;'
         type *, '  Enter 1 to change sign of y axis in input file:'
         accept *, iysign
      endif         

      type *, 'Enter 1 if input file scale is not in arcsec:'
      accept *, iscale
      if (iscale .eq. 1) then
         type *, 'Enter pixel scale (in arcsec/pixel):'
         accept *, platescale
      endif

      type *, 'Enter number of lines to skip at top of the input file:'
      accept *, iskip

c
c open the file and read its contents
c
c the flag is = 0 for data to be used; =1 for data to be excluded
c

      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,k
         read(11,*,end=1001,err=1001) inumobj(i),xobj(i),
     +        yobj(i),wslit(i),magin(i)
         inumin=inumin+1
         iflag(i)=0
         if (isense .eq. 1) then
            if (ixsign .eq. 1) then
               xobj(i)=-1. * xobj(i)
            endif
            if (iysign .eq. 1) then
               yobj(i)=-1. * 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.
      enddo
 1001 close(unit=11,status='keep')
      inumnow=inumin

c
c give opportunity to change the general size of the plot and the
c     enlargement factor
c

      type *, 'Enter 1 to accept usual parameters for size of plot'
      type *, '  and enlargement factor (720"x480";mag=4);'
      type *, '  2 to change:'
 199  accept *, isize
      if (isize .eq. 2) then
 198     type *, 'Enter xsize/"(real), ysize/"(real), magfactor(real):'
         accept *, xsize, ysize, enlarge

         if ( xsize/scale*enlarge/25.4 .gt. 10.5) then
            type *, 'The plot is too large for the page; choose',
     +           ' a new xsize'
            goto 198
         endif
         if ( xsize/scale*enlarge/25.4 .gt. 8.) then
            type *, 'The plot is too large for the page; choose',
     +           ' a new ysize'
            goto 198
         endif

      else if (isize .eq. 1) then
         xsize=720.
         ysize=480.
         enlarge=4.
      else 
         type *, 'Enter 1 or 2:'
         goto 199
      endif

c
c this is the general loop for reading unflagged data from the input
c     list into the working vectors
c

 2000 icount=0
      x0=0.
      y0=0.

      do i=1,inumin

         if (iflag(i) .eq. 0) then
            icount=icount+1
            xproj(icount)=xobj(i)*cos(theta) + yobj(i)*sin(theta)
            yproj(icount)=yobj(i)*cos(theta) - xobj(i)*sin(theta)
            wproj(icount)=wslit(i)
            iobjproj(icount)=inumobj(i)
            slitposeproj(icount)=slitpose(i)
            slitposwproj(icount)=slitposw(i)
            magproj(icount)=magin(i)

            x0 = x0 + xproj(icount)
            y0 = y0 + yproj(icount)

         endif

      enddo

c
c put object numbers of zero into the unused ones
c

      do i=icount+1,inumin
         iobjproj(i)=0
      enddo

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

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

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

c      type 110, x0, y0
      type 110, x0unproj, y0unproj

c
c sort the objects into ascending x order, while keeping the other arrays
c     in the associated positions
c piksr7 subroutine is *adapted* from numerical recipes; 4 indicates that
c     there are six dependent vectors (yproj,wproj,iobjproj,magproj,
c     slitposeproj,slitposwproj) on the one sorting vector (xproj)
c

      call piksr7(icount,xproj,yproj,wproj,iobjproj,slitposeproj,
     +     slitposwproj,magproj)

c
c make all (x,y) positions relative to the centroid
c

      do i=1,icount
         xproj(i) = xproj(i) - x0 + xoff
         yproj(i) = yproj(i) - y0 + yoff
         if (xproj(i) .lt. (xsize/(-2.)-15.) 
     +        .or. xproj(i) .gt. (xsize/2.+15.)
     +        .or. yproj(i) .lt. (ysize/(-2.)-15.)
     +        .or. yproj(i) .gt. (ysize/2.+15.) ) then
            type 130, iobjproj(i), xproj(i),yproj(i)
         endif
         if (yproj(i) .gt. 60. .or. yproj(i) .lt. -60.) then
            type 131, iobjproj(i), xproj(i),yproj(i)
         endif

      enddo

c
c set up the slit coordinates
c

      do i=1,icount
         if (i .eq. 1) then
            xstart(i) = xproj(i) - 15.
         else
c remove one pixel = 0.4" from length on each end of slit
            xstart(i)=(xproj(i)+xproj(i-1))/2.+slitposeproj(i)+0.4
         endif
         if (i .eq. icount) then
            xend(i) = xproj(i) + 15.
         else
            xend(i) = (xproj(i)+xproj(i+1))/2.+slitposwproj(i)-0.4
         endif
         if (xproj(i)-xstart(i) .lt. 3.) then
            type 170, iobjproj(i)
         endif
         if (xend(i)-xproj(i) .lt. 3.) then
            type 171, iobjproj(i)
         endif
c
c modification:  now gives the percentage of the slit for the object
c     measured from the E to the W (i.e. left to right on the plot)
c
c         percent(i) = (xproj(i)-xstart(i)) / (xend(i)-xstart(i))
         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)

         if (lproj(i) .le. 10.) then
            type 172, iobjproj(i)
         endif

         ystart(i) = yproj(i) - wslit(i)/2.
         yend(i) = yproj(i) + wslit(i)/2.

      enddo
c
c plot the current working set of vectors
c

      if (ioption .ne. 0) then
         call pgend
      endif

      if (ioption .eq. 6) then
c /ps to create pgplot.ps 
         call pgbegin(0,'/ps',1,1)
      else 
c /xterm to terminal
         type *, 'Enter /xw for xwindow; /xt for xterm:'
         call pgbegin(0,'?',1,1)
      endif
c
c size of plot in inches; x is 720 arcsec, y is 480 arcsec
c
      x1=0.25
      x2=xsize/scale*enlarge/25.4+0.25
      y1=0.25
      y2=ysize/scale*enlarge/25.4+0.25
      call pgvsize(x1,x2,y1,y2)
c
c note:  x +ve in the right direction = E (sky parity)
c     y +ve in the upper direction = N
c

      x1=xsize/(-2.)
      x2=xsize/(2.)
      y1=ysize/(-2.)
      y2=ysize/(2.)
c
c modification:  now N up and E left
c
c      call pgwindow(x1,x2,y1,y2)
      call pgwindow(x2,x1,y1,y2)
      call pgslw(1)
      call pgbox('BC',0.,0.,'BC',0.,0.)

c
c enter the name in the upper LH corner
c

      call pgsch(0.6)
      call pgtext(x2-20.,y2-12.,string1)

      write(string2,150) enlarge
      call pgtext(x2-20.,y2-24.,string2)

      theta=theta*360./2./pi
c      write(string3,151) theta
      write(string3,155) theta+180.
      theta=theta/360.*2.*pi
      call pgtext(x2-20.,y2-36.,string3)
      write(string4,152) x0unproj,y0unproj

      call pgtext(x2-20.,y2-48.,string4)
      call pgsch(1.0)

c
c plot the rectangular slits 
c

      do i=1,icount
         xpts(1)=xstart(i)
         ypts(1)=ystart(i)
         xpts(2)=xend(i)
         ypts(2)=ystart(i)
         xpts(3)=xend(i)
         ypts(3)=yend(i)
         xpts(4)=xstart(i)
         ypts(4)=yend(i)
         xpts(5)=xstart(i)
         ypts(5)=ystart(i)
c         call pgslw(3)
c number of line strokes
         call pgslw(1)

c
c psgsfs(1=solid;2=hollow)
c

         call pgsfs(1)

         call pgpoly(5,xpts,ypts)
         call pgslw(1)
      enddo

c
c if plotting to the terminal, plot the number of each slit above it
c
      
      if (ioption .ne. 6 .or. inumbers .eq. 1) then
         call pgsch(0.8)
         do i=1,icount
            text='    '
            write(text,140) iobjproj(i)
            x1=xproj(i)
            y1=yproj(i)+10.
            if (iobjproj(i) .ge. 1000) then
               call pgptext(x1,y1,0.,0.5,text)
            else if (iobjproj(i) .ge. 100) then
               call pgptext(x1,y1,0.,0.5,text(2:4))
            else if (iobjproj(i) .ge. 10) then
               call pgptext(x1,y1,0.,0.5,text(3:4))
            else 
               call pgptext(x1,y1,0.,0.5,text(4:4))
            endif
            
            y1=yproj(i)
            call pgpoint(1,x1,y1,2)

         enddo
         call pgsch(1.0)
      endif
      inumbers=2

      iprevious = ioption

c
c necessary for the pgplot.ps output to end the plot
c

      if (ioption .eq. 6) then
         call pgend
      endif

c
c query the user if they desire a change in any of the parameters
c

 4000 type *, 'Enter 1 to change the rotation angle,'
      type *, '      2 to exclude an object,'
      type *, '      3 to include an object,'
      type *, '      4 to plot again on the terminal,'
      type *, '      5 to type the object list to the terminal,'
      type *, '      6 to write an output file,'
      type *, '      7 to edit lengths of slits'
      type *, '      8 to include all objects'
      type *, '      9 to exclude all objects'
      type *, '      10 to add overall offset in x or y to all slits:'
      type *, '      11 to exit program, no output:'
      type *, ' '
      accept *, ioption

      if (ioption .eq. 1) then
         type *, 'Enter a rotation angle (CCW from input x,y; degrees):'
         accept *, theta
         theta=theta/360.*2.*pi
         type *, 'Note that re-projection onto this theta will not be'
         type *, '  calculated until you replot the data...'
         goto 4000
      else if (ioption .eq. 2) then
         type *, 'Enter number of object to exclude:'
         accept *, iex
         do i=1,inumin
            if (inumobj(i) .eq. iex) then
               iflag(i)=1
               goto 4001
            endif
         enddo
         type *, 'No object found with that number, try again:'
 4001    goto 4000
      else if (ioption .eq. 3) then
         type *, 'Enter number of object to include:'
         accept *, iinc
         do i=1,inumin
            if (inumobj(i) .eq. iinc) then
               iflag(i)=0
               goto 4002
            endif
         enddo
         type *, 'No object found with that number, try again:'
 4002    goto 4000
      else if (ioption .eq. 4) then
         goto 2000
      else if (ioption .eq. 5) then
         type 103, inumnow
         type *, 'Included objects:'
         type *, 'Object   x/"    y/"   width/"    length/"    percent',
     +        '  L/"(Wside)  L/"(Eside)'
         type *, '                                             from top'
         type *, '  '
         do i=1,inumin
            if (iflag(i) .eq. 0) then
               do j=1,inumnow
                  if (inumobj(i) .eq. iobjproj(j)) then
                     type 102, iobjproj(j),xproj(j),yproj(j),wproj(j),
     +                    lproj(j),percent(j),lwest(j),least(j)
                  endif
               enddo
            endif
         enddo
         type *, '  '
         type *, 'Excluded objects:'
         type *, 'Object   x/"    y/"   width/"'
         type *, '  '
         do i=1,inumin
            if (iflag(i) .eq. 1) then
               type 104, inumobj(i),xobj(i),yobj(i),wslit(i)
            endif
         enddo
         goto 4000
      else if (ioption .eq. 6) then
         type *, 'Enter output filename for textfile output:'
         accept *, outfile
         open(unit=12,file=outfile,form='formatted',status='new')
         theta=theta*360./2./pi
c
c modification:  more notation at top of file
c
c         write(12,180) x0unproj,y0unproj,theta,180.+theta
c         type 180, x0unproj,y0unproj,theta,180.+theta

         write(12,181) x0unproj,y0unproj
         write(12,182) theta,theta+180.
         write(12,183)
         type 181,x0unproj,y0unproj
         type 182,theta,theta+180.
         theta=theta/360.*2.*pi

         do i=1,inumin
            if (iflag(i) .eq. 0) then
               do j=1,inumnow
                  if (inumobj(i) .eq. iobjproj(j)) then
c            write(12,101) iobjproj(j),xproj(j)+x0,yproj(j)+y0,wproj(j),
c     +           lproj(j),percent(j),xproj(j),yproj(j),least(j),lwest(j)
                     write(12,101) iobjproj(j),xobj(i),yobj(i),
     +                    wproj(j),magproj(j),lproj(j),percent(j),
     +                    xproj(j),yproj(j),least(j),lwest(j)
                  endif
               enddo
            endif
         enddo
         type *, 'Postscript version of plot written to pgplot.ps'
         type *, 'Print by "lpr -s -Pps pgplot.ps"'
         type *, ' '
         close(unit=12,status='keep')
         type *, 'Enter 1 to include object numbers in plot'
         type *, '    (for reference purposes),'
         type *, '  2 for no numbers (for the actual mask):'
         accept *, inumbers
         goto 2000
      else if (ioption .eq. 7) then
         type *, 'Enter object number on W side of slit interface:'
         accept *, iobjw
         type *, 'Enter object number on E side of slit interface:'
         accept *, iobje
         type *, 'Enter amount to change interface (E +ve):'
         accept *, 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
         goto 4000
      else if (ioption .eq. 8) then
         do i=1,inumin
c flag=0 include; flag=1 exclude
            iflag(i)=0
         enddo
         goto 4000
      else if (ioption .eq. 9) then
         do i=1,inumin
            iflag(i)=1
         enddo
         goto 4000
      else if (ioption .eq. 11) then
         goto 5000
      else if (ioption .eq. 10) then
         type *, 'Enter total offset from centroid in x,y (arcsec)'
         accept *, xoff,yoff
         goto 4000
      else
         type *, 'Enter value 1-11 only:'
         goto 4000
      endif

 5000 continue

      call pgend

      stop

 101  format(i4,x,2(f7.2,x),2(f5.2,x),f7.2,x,f5.3,x,4(f7.2,x))
 102  format(i4,x,2(f7.2,x),f5.2,x,f7.2,x,f5.3,x,2(f7.2,x))
 103  format('Number of objects in current plot:',3x,i4)
 104  format(i4,x,2(f7.2,x),f5.2)
 110  format('Centroid of slits is at x=',f7.2,' arcsec,  y=',
     +     f7.2,' arcsec')
 120  format('Position angle = ',f7.2,' degrees, centroid x0=',f7.2,
     +     ' arcsec, y0=',f7.2)
 130  format('Object number ',i4,' at (',f7.2,x,f7.2,
     +     ') is out of the field of view')
 131  format('Object number ',i4,' at (',f7.2,x,f7.2,
     +     ') is out of the recommended range of |y|<60"')
 140  format(i4)
 150  format('Magnification=',f5.2)
 151  format('E=',f7.2,' deg CW from +ve x')
 152  format('x0=',f7.2,'" y0=',f7.2,'" relative to input (x=0,y=0)')
 155  format('COSMIC position angle = ',f7.2,' degrees')
 160  format('Theta=',f5.2,' deg CCW from x +ve=E')
 170  format('Object #',i4,' has a slit on W side of object < 3"')
 171  format('Object #',i4,' has a slit on E side of object < 3"')
 172  format('Object #',i4,' has a slit < 10"')
 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)
 181  format('Centroid x0=',f8.2,'" and y0=',f8.2,'" relative to the',
     +     ' input (x=0,y=0) coordinates')
 182  format('Rotation angle from input coordinates:',f8.2,'  COSMIC ',
     +     'base position angle = ',f8.2)
 183  format('Obj#   Xobj/"    Yobj/"    Slitwidth/"   Mag ',
     +     '  Slitlength/"  %(along slit E->W)  Xrel/"   Yrel/"',
     +     '  SlitlengthW/"   SlitlengthE/"')

      end



      subroutine piksr7(n,arr,b1rr,b2rr,b3rr,b4rr,b5rr,b6rr)
      dimension arr(n),b1rr(n),b2rr(n),b3rr(n),b4rr(n),b5rr(n),b6rr(n)
      do j=2,n
         a=arr(j)
         b1=b1rr(j)
         b2=b2rr(j)
         b3=b3rr(j)
         b4=b4rr(j)
         b5=b5rr(j)
         b6=b6rr(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)
            b4rr(i+1)=b4rr(i)
            b5rr(i+1)=b5rr(i)
            b6rr(i+1)=b6rr(i)
         enddo
         i=0
 10      arr(i+1)=a
         b1rr(i+1)=b1
         b2rr(i+1)=b2
         b3rr(i+1)=b3
         b4rr(i+1)=b4
         b5rr(i+1)=b5
         b6rr(i+1)=b6
      enddo
      return
      end