c***************************************** program moemask c c maintained and partially written by A.McWilliam 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 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 Mask PA 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 "moe.common.blocks" when compiling. c The file contains the following statements: c c--------------start of moe.common.blocks--------------------------------- c c integer kmax c parameter (kmax=4000) ! maximum number of objects in a mask list c common/iindata/ inumin,inumobj c integer inumin,inumobj(kmax) c common/islit/ iobjproj,iflag,kcount,icount,iscale c integer iobjproj(kmax),iflag(kmax),kcount(kmax),icount,iscale c common/iparams/ isense,ixsign,iysign,iwarn,ioffset,kpri,iwfccd, c % iskip c integer isense,ixsign,iysign,iwarn,ioffset,kpri,iwfccd,iskip c common/rindata/xobj,yobj,priority,xmag,pa,wslit,slitposw,slitpose, c . orah,oram,oras,odecd,odecm,odecs,ora,odec,dra_as,ddec_as c real*8 xobj(kmax),yobj(kmax),priority(kmax),xmag(kmax),pa(kmax) c real*8 wslit(kmax),slitposw(kmax),slitpose(kmax) c c Note: ora odec etc are the object ra and dec values c c real*8 orah(kmax),oram(kmax),oras(kmax),odecd(kmax),odecm(kmax) c real*8 odecs(kmax),ora(kmax),odec(kmax),dra_as(kmax),ddec_as(kmax) c common/cindata/ora_str,odec_str c character*12 ora_str(kmax),odec_str(kmax) c common/rslit/ xproj,yproj,xstart,xend,ystart,yend,slitposwproj, c % slitposeproj,x0,y0,x0unproj,y0unproj,wproj,lproj, c % paproj,least,lwest,percent,projpriority,xspecr_ccd, c % xspecb_ccd,xmask_ccd,ymask_ccd,xmm,ymm,xmmproj,ymmproj c real*8 xproj(kmax),yproj(kmax),xstart(kmax),xend(kmax),ystart(kmax) c real*8 yend(kmax),slitposwproj(kmax),slitposeproj(kmax),x0,y0 c real*8 x0unproj,y0unproj,wproj(kmax),lproj(kmax),least(kmax) c real*8 paproj(kmax),lwest(kmax),percent(kmax),projpriority(kmax) c real*8 xspecr_ccd(kmax),xspecb_ccd(kmax),xmask_ccd(kmax),ymask_ccd(kmax) c real*8 xmm(kmax),ymm(kmax),xmmproj(kmax),ymmproj(kmax) c common/rslit2/ shiftx,shifty c real*8 shiftx,shifty c common/rparams/ fl,asperunit,scale,wfpix,starwidth,xsize,ysize,ystack, c % xccd_size,endpix,endslit,slitmin,enlarge,xoff,yoff,globslit, c % holediam,thinrot,theta,thetad,maskpa,slitpa,cosmic_scale, c % red_limit,blue_limit c real*8 fl,asperunit,scale,wfpix,starwidth,xsize,ysize,ystack,xccd_size,endpix c real*8 endslit,slitmin,enlarge,xoff,yoff,globslit,holediam,thinrot,theta c real*8 thetad,maskpa,slitpa,cosmic_scale,red_limit,blue_limit c common/lparams/auto_priority,auto_offset,tilted_slits,holes,stacked, c % cosmic,wfccd,interactive,sw074,sw103,sw125,sw147, c % sw177,sw199,nod_shuffle c logical auto_priority, auto_offset, tilted_slits, holes, stacked, c % cosmic, wfccd, interactive, sw074, sw103, sw125, sw147, c % sw177, sw199, nod_shuffle c common/cparams/ title,infile,psfile,gmoutput,exclfile,stdout c character*64 title*32,infile,psfile,gmoutput,exclfile,stdout c c Note: the following coordinates are for the field center c c common/rcoords/ rah,ram,ras,decd,decm,decs,epoch,rahproj,ramproj, c % rasproj,decdproj,decmproj,decsproj,dec_sign,ra_as,dec_as c real*8 rah,ram,ras,decd,decm,decs,epoch,rahproj,ramproj,rasproj c real*8 decdproj,decmproj,decsproj,dec_sign,ra_as,dec_as c common/lcoords/coords c logical coords c *** DEBUG c common/dbg/debug c logical debug c c c--------------end of moe.common.blocks----------------------------------- c implicit real*8 (a-h,o-z) include 'moe.common.blocks' c ** DEBUG data debug/.false./ c c data scale,wfpix/2.902,0.380/ ! original telescope scale computed by B.Sutin c data scale,wfpix/2.880,0.3767/ ! Baade scale measured my J.Mulchaey data fl,scale,wfpix/71161.1,2.89856,0.1110/ ! Baade IMACS fl, mask scale, ccd scale data starwidth/5.80/ ! width of the alignment star boxes in arcsec data xsize,ysize,ystack/926.0,1200.0,900.0/ !limits of the imaged mask field in arcsec c data xsize,ysize,ystack/926.0,600.0,900.0/ !limits of the imaged mask field in arcsec data xccd_size/8342.0/ data endpix,endslit,slitmin,globslit,holediam/500.0,22.5,5.0,0.50,1.94/ data enlarge/1.05/ ! enlargement factor data xoff,yoff,shiftx,shifty/0.0,0.0,0.0,0.0/ ! intitialize some parameters... data holes,auto_priority,tilted_slits,coords/.false.,.true.,.false.,.false./ data auto_offset/.false./ data stacked,cosmic,wfccd,interactive/.false.,.false.,.true.,.true./ data sw074, sw103, sw125, sw147, sw177, sw199 /6*.false./ data isense,iscale,iwfccd,iskip,ioffset,kpri/0,0,0,0,0,0/ call read_par_file c call query_user ! enter input file name, title and set parameters call read_data call convert_coords call select_all_objects(iflag,inumin) !initialize all objects to be considered call create_slit_mask ! select objects and create slit mask information c if (interactive.eq..true.) then if (interactive) then call make_xwindow_spectrum_plot ! plot slits and spectrm on screen call interact_with_user call pgend endif call produce_output_files stop ' NORMAL STOP ' end character*(64) function literal(string) character*(*) string nchar = len(string) do i = 1,nchar if (string(i:i).ne.' ') then istart = i go to 100 endif enddo 100 continue do j = nchar,1,-1 if (string(j:j).ne.' ' ) then iend = j go to 200 endif enddo 200 continue literal = string(istart:iend) return end subroutine read_par_file implicit real*8 (a-h,o-z) include 'moe.common.blocks' character*102 line character*(64) literal logical arc_scale parameter (raddeg= 57.29577951) open(unit=23,file='moe.par',status='old') do i = 1,100 read(23,'(a)',end=900)line if (line(:21).eq.'input ') then c read(line(23:),'(a)')infile infile = literal(line(23:)) elseif (line(:21).eq.'gm output ') then gmoutput = literal(line(23:)) elseif (line(:21).eq.'postscript output ') then psfile = literal(line(23:)) elseif (line(:21).eq.'unused star file name') then exclfile = literal(line(23:)) if (exclfile.ne.' ') % open(unit=40,file=exclfile,status='unknown') elseif (line(:21).eq.'standard output file ') then stdout = literal(line(23:)) if (stdout.ne.' ') % open(unit=41,file=stdout,status='unknown') elseif (line(:21).eq.'plot title ') then read(line(23:),'(a32)')title elseif (line(:21).eq.'RA ') then call get_ra(line(23:)) elseif (line(:21).eq.'DEC ') then call get_dec(line(23:)) elseif (line(:21).eq.'Epoch ') then read(line(23:),*)epoch elseif (line(:21).eq.'Cosmic ') then read(line(23:),*)cosmic elseif (line(:21).eq.'Alignment size ') then read(line(23:),*)starwidth elseif (line(:21).eq.'Positions in arcsec? ') then read(line(23:),*)arc_scale if (.not.arc_scale)iscale = 1 elseif (line(:21).eq.'Arc sec per unit ') then read(line(23:),*)asperunit elseif (line(:21).eq.'Mask scale (as/mm)') then read(line(23:),*)scale elseif (line(:21).eq.'Pixel scale (as/pix)') then read(line(23:),*)wfpix elseif (line(:21).eq.'skip lines ') then read(line(23:),*)iskip elseif (line(:21).eq.'Use auto Offsets? ') then read(line(23:),*) auto_offset elseif (line(:21).eq.'X offset (arcsec) ') then read(line(23:),*)xoff elseif (line(:21).eq.'Y offset (arcsec) ') then read(line(23:),*)yoff elseif (line(:21).eq.'Auto Priority ') then read(line(23:),*) auto_priority elseif (line(:21).eq.'Field X-size (arcsec)') then read(line(23:),*)xsize elseif (line(:21).eq.'Field Y-size (arcsec)') then read(line(23:),*)ysize elseif (line(:21).eq.'Holes ') then read(line(23:),*)holes c elseif (line(:21).eq.'Hole diameter ') then c read(line(23:),*)holediam c elseif (line(:21).eq.'Slits? ') then c read(line(23:),*)slits elseif (line(:21).eq.'Interactive? ') then read(line(23:),*)interactive elseif (line(:21).eq.'Slit width (as) ') then read(line(23:),*) globslit elseif (line(:21).eq.'Slit width 0.74 as ') then read(line(23:),*) sw074 if (sw074) globslit = 0.74 elseif (line(:21).eq.'Slit width 1.03 as ') then read(line(23:),*) sw103 if (sw103) globslit = 1.03 elseif (line(:21).eq.'Slit width 1.25 as ') then read(line(23:),*) sw125 if (sw125) globslit = 1.25 elseif (line(:21).eq.'Slit width 1.47 as ') then read(line(23:),*) sw147 if (sw147) globslit = 1.47 elseif (line(:21).eq.'Slit width 1.99 as ') then read(line(23:),*) sw199 if (sw199) globslit = 1.99 elseif (line(:21).eq.'min slit sep. (pix) ') then read(line(23:),*)endpix elseif (line(:21).eq.'min end-object (as) ') then read(line(23:),*)endslit elseif (line(:21).eq.'min slit length (as) ') then read(line(23:),*)slitmin elseif (line(:21).eq.'Tilted slits? ') then read(line(23:),*)tilted_slits elseif (line(:21).eq.'Nod and Shuffle dirn?') then read(line(23:),*)nod_shuffle if (nod_shuffle) then xccd_width = 8242.0 else xccd_width = 8342.0 endif elseif (line(:21).eq.'Stack spectra? ') then read(line(23:),*)stacked elseif (line(:21).eq.'min stack dist. (pix)') then read(line(23:),*)ystack elseif (line(:21).eq.'Red Limit (A) ') then read(line(23:),*)red_limit elseif (line(:21).eq.'Blue Limit (A) ') then read(line(23:),*)blue_limit elseif (line(:21).eq.'Slit PA ') then read(line(23:),*) slitpa if (slitpa.ge. 0.000 .and. slitpa.le.360.0) then maskpa = slitpa - 90.0 thinrot = maskpa thetad = thinrot if (maskpa.lt.0.0) maskpa = maskpa + 360.0 theta = thetad/raddeg else write(6,'(a)')'Error: Slit PA must be between 0 and 360' write(6,'(a)')line stop endif else write(6,'(a)')'Error: unrecognized parameter' write(6,'(a)')line stop endif enddo 900 continue if (stacked) then ystack = ystack*wfpix else ystack = ysize endif if (holes) then c globslit = holediam endslit = globslit/2.0 slitmin = globslit endif call open_gm_files(globslit,gmoutput) return end subroutine open_gm_files(sw,gmoutput) c c decide on single or double pass and open the GM code output files c implicit real*8 (a-h,o-z) real*8 sw character*64 gmoutput,gmout1*70,gmout2*70 do j=1,64 if (gmoutput(j:j).eq.' ') goto 1 enddo 1 len=min(j-1,8) if (sw.lt.0.90) then gmout1 = gmoutput(1:len) // '.ld1' open (unit=30,file=gmout1,status='unknown') gmout2 = gmoutput(1:len) // '.ld2' open (unit=20,file=gmout2,status='unknown') else gmout2 = gmoutput(1:len) // '.lds' open (unit=20,file=gmout2,status='unknown') endif return end subroutine interact_with_user c c allow user to modify parameters and select options interactively c implicit real*8 (a-h,o-z) include 'moe.common.blocks' logical numbers integer ioption numbers = .true. do while(.true.) call get_commands(ioption) if (ioption .eq. 1) then call change_slit_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 c call prompt_for_numbers(numbers) c dxpaper = 16.5 c dypaper = 10.5 call make_postscript_file call rename_photomask_pgplot_file elseif (ioption .eq. 7) then call make_postscript_spectrum_file call rename_photomask_pgplot_file 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 call make_xwindow_plot elseif (ioption .eq. 12) then return else write(6, '(a)')'Enter value 1-12 only:' call get_commands(ioption) endif if (ioption.le.6 .or. ioption.eq.10) then call create_slit_mask call make_xwindow_spectrum_plot endif enddo return end c subroutine multisort(n,arr,arb,axrr,ayrr,b0rr,b1rr,b2rr,b3rr,nb4rr,b5rr, c % b6rr,nb7rr,b8rr) subroutine multisort(n,arr,arb,axrr,ayrr,axm,aym,b0rr,b1rr,b2rr,b3rr,nb4rr,b5rr, % b6rr,nb7rr,b8rr) implicit real*8 (a-h,o-z) 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*8 arr(kmax),arb(kmax),b0rr(kmax),b1rr(kmax),b2rr(kmax),b3rr(kmax) real*8 b5rr(kmax),b6rr(kmax),b8rr(kmax),axrr(kmax),ayrr(kmax) real*8 axm(kmax),aym(kmax) integer nb4rr(kmax),nb7rr(kmax) do j=2,n a=arr(j) ab=arb(j) ax=axrr(j) ay=ayrr(j) axmm=axm(j) aymm=aym(j) b0=b0rr(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) arb(i+1)=arb(i) axrr(i+1)=axrr(i) ayrr(i+1)=ayrr(i) axm(i+1)=axm(i) aym(i+1)=aym(i) b0rr(i+1)=b0rr(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 arb(i+1)=ab axrr(i+1)=ax ayrr(i+1)=ay axm(i+1)=axmm aym(i+1)=aymm b0rr(i+1)=b0 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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' c 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)*dcos(theta) + yobj(iadd)*dsin(theta) + shiftx yproj(icount)=yobj(iadd)*dcos(theta) - xobj(iadd)*dsin(theta) + shifty c c check for out of bounds slit c if (xproj(icount).gt.0.5*xsize .or. xproj(icount).lt.-0.5*xsize .or. % yproj(icount).gt.0.5*ysize .or. yproj(icount).lt.-0.5*ysize ) then write(6,'(15hError: object #,i10,22h is out of mask bounds)')iinc iflag(iadd) = 1 icount = icount - 1 return endif 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 = j xlow = xdist endif else if (xdist.le.xup) then iup = j xup = xdist endif endif enddo c c account for adjacent align boxes c if (projpriority(ilow).lt.0.0) xlow = xlow + starwidth/2.0 if (projpriority(iup) .lt.0.0) xup = xup - starwidth/2.0 c c now find conflicts with adjacent slits and suggest deletions c c stars too close xstart_up = (xproj(icount)+xproj(iup))*0.5 + wfpix*endpix*0.5 xend_low = (xproj(icount)+xproj(ilow))*0.5 - wfpix*endpix*0.5 slen_up = xend(iup) - xstart_up slen_low = xend_low - xstart(ilow) dmin = endpix*wfpix + 2.0*endslit slength = xup - xlow - dmin if (-xlow.lt.dmin .or. slen_low.lt.slitmin) then write(6,'(22h Conflict with object ,i10)')iobjproj(ilow) call delete(iobjproj(ilow),iflag(kcount(ilow))) endif if (xup.lt.dmin .or. slen_up.lt.slitmin) then write(6,'(22h Conflict with object ,i10)')iobjproj(iup) call delete(iobjproj(iup),iflag(kcount(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)')iobjproj(ilow),iobjproj(iup) call delete(iobjproj(ilow),iflag(kcount(ilow))) call delete(iobjproj(iup),iflag(kcount(iup))) endif return end integer function found_object(iinc) implicit real*8 (a-h,o-z) include 'moe.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 implicit real*8 (a-h,o-z) include 'moe.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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' call get_input_file_name(infile) call get_title(title) call get_coords call set_alignment_warning call change_defaults call change_slit_pa call set_global_slitwidth return end subroutine read_data implicit real*8 (a-h,o-z) include 'moe.common.blocks' character*80 line integer left(6),right(6),nstr 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 n = 0 do i = 1,kmax n = n + 1 read(11,'(a)',end=1001, err=1001) line call parse_string(line,left,right,nstr) c program star if (line(1:1).eq.'@' .and. nstr.ge.5) then read(line(left(2):right(2)),*)inumobj(n) ora_str(n)=line(left(3):right(3)) odec_str(n)=line(left(4):right(4)) read(line(left(5):right(5)),*)xmag(n) priority(n) = xmag(n) wslit(n) = globslit if (nstr.ge.6) read(line(left(6):right(6)),*)pa(n) c alignment star elseif (line(1:1).eq.'*' .and. nstr.ge.5) then read(line(left(2):right(2)),*)inumobj(n) ora_str(n)=line(left(3):right(3)) odec_str(n)=line(left(4):right(4)) read(line(left(5):right(5)),*)xmag(n) priority(n) = -4.0 wslit(n) = starwidth c comment line elseif (line(1:1).eq.'!') then n = n - 1 inumin = inumin - 1 goto 900 endif inumin = inumin + 1 if (inumin.ge.kmax) stop ' too many entries in input list' slitpose(n)=0. slitposw(n)=0. 900 continue enddo 1001 close(unit=11,status='keep') return end subroutine parse_string(line,left,right,nstr) c a routine to count the number of sub-strings and their locations within a line implicit real*8 (a-h,o-z) character*80 line integer left(6),right(6),nstr logical in_word nstr = 0 in_word = .false. do i = 1,80 if (in_word) then if (line(i:i).eq.' ') then right(nstr) = i-1 in_word = .false. endif else if (line(i:i).ne.' ') then in_word=.true. nstr = nstr + 1 left(nstr) = i endif endif enddo return end c subroutine read_data c c implicit real*8 (a-h,o-z) c include 'moe.common.blocks' c c open(unit=11,file=infile,form='formatted',status='old') c inumin=0 c c if (iskip .gt. 0) then c do i=1,iskip c read(11,*,end=1001,err=1001) c enddo c endif c c do I=1,kmax c if (tilted_slits) then c read(11,*,end=1001,err=1001) inumobj(i),xobj(i), c % yobj(i),priority(i),xmag(i),pa(i) c else c read(11,*,end=1001,err=1001) inumobj(i),xobj(i), c % yobj(i),priority(i),xmag(i) c endif c wslit(i) = globslit c inumin=inumin+1 c if(inumin.ge.kmax) stop ' too many entries in input list' c c if (isense .eq. 1) then c if (ixsign .eq. 1) then c xobj(i) = -xobj(i) c endif c if (iysign .eq. 1) then c yobj(i) = -yobj(i) c endif c endif c if (iscale .eq. 1) then c xobj(i)=xobj(i)*asperunit c yobj(i)=yobj(i)*asperunit c endif c slitpose(i)=0. c slitposw(i)=0. cc ... if this object is a star, set its slitwidth to "starwidth' c if(priority(i).lt.0.0) wslit(i) = starwidth c enddo c1001 close(unit=11,status='keep') c c return c end subroutine convert_coords c convert from ra/dec strings to numbers for the object coordinates implicit real*8 (a-h,o-z) include 'moe.common.blocks' parameter (raddeg= 57.29577951) integer i1,i2,i3 do i = 1,inumin i1 = index(ora_str(i),':') read(ora_str(i)(:i1-1),*)orah(i) i2 = index(ora_str(i)(i1+1:),':') read(ora_str(i)(i1+1:i1+i2-1),*)oram(i) i3 = index(ora_str(i)(i1+i2+1:),':') read(ora_str(i)(i1+i2+i3+1:),*)oras(i) i1 = index(odec_str(i),':') read(odec_str(i)(:i1-1),*)odecd(i) i2 = index(odec_str(i)(i1+1:),':') read(odec_str(i)(i1+1:i1+i2-1),*)odecm(i) i3 = index(odec_str(i)(i1+i2+1:),':') read(odec_str(i)(i1+i2+i3+1:),*)odecs(i) c now compute RA and DEC in arcsec and then get difference from field center c watch-out for fields including a pole (N or S)! An error will occur ora_as = nint(orah(i)) * 3600.0 + nint(oram(i))*60.0 + oras(i) odec_as = abs(nint(odecd(i)))*3600.0 + nint(odecm(i))*60.0 + odecs(i) osign = nint(abs(odecd(i))/odecd(i)) odec_as = osign * odec_as odec(i) = odec_as/3600.0d0 ora(i) = ora_as/3600.0d0 dra = ora_as - ra_as if (int(rah).eq.23) then if (int(orah(i)).eq.0) then dra = dra + 86400.0 endif elseif (int(rah).eq.0) then if (int(orah(i)).eq.23) then dra = dra - 86400.0 endif endif dra_as(i) = 15.0 * dra * dcos((odec_as/3600.0)/raddeg) c if (dra_as(i).lt.-36000.0) then c dra_as(i) = dra_as(i) + 1296000.0 c elseif (dra_as(i).gt.36000.0) then c dra_as(i) = dra_as(i) - 1296000.0 c endif ddec_as(i) = odec_as - dec_as c xobj(i) = dra_as(i) xobj(i) = -dra_as(i) ! -ve because E is to left in slit plane yobj(i) = ddec_as(i) xo = xobj(i) yo = yobj(i) call compute_slit_xymm(xo,yo,odec(i),dec_as,fl,xmm(i),ymm(i)) enddo return end subroutine make_new_xylist c c compute projected x,y positions of un-deleted objects and centroid c implicit real*8 (a-h,o-z) include 'moe.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 xo = xmm(i)*scale yo = ymm(i)*scale call project_to_pa(xo,yo,theta,xproj(icount),yproj(icount)) xmmproj(icount) = -xproj(icount)/scale ! -ve for optics reversal? ymmproj(icount) = -yproj(icount)/scale ! -ve for optics reversal? 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*dcos(theta) - y0*dsin(theta) y0unproj=y0*dcos(theta) + x0*dsin(theta) c write(6,'(a)')'x0,y0' c write(6,*)x0,y0 c write(6,'(a)')'x0unproj,y0unproj,theta' c write(6,*) x0unproj,y0unproj,theta return end subroutine project_to_pa(xo,yo,theta,xp,yp) c c computes the projected x,y mask positions for a given position angle c implicit real*8 (a-h,o-z) parameter (raddeg= 57.29577951d0) c c apply the transform c xp = dcos(theta)*xo + dsin(theta)*yo yp = dcos(theta)*yo - dsin(theta)*xo return end subroutine sort_slits c c sort slit information by ascending xproj c implicit real*8 (a-h,o-z) include 'moe.common.blocks' call multisort(icount,xspecr_ccd,xspecb_ccd,xmask_ccd,ymask_ccd, % xmmproj,ymmproj,xproj,yproj,paproj,wproj,iobjproj, % slitposeproj,slitposwproj,kcount,projpriority) return end subroutine shift_xy implicit real*8 (a-h,o-z) include 'moe.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 (auto_offset) then shiftx = -(xoff*dcos(theta)+yoff*dsin(theta)) - x0 shifty = -(yoff*dcos(theta)-xoff*dsin(theta)) - y0 else shiftx = -(xoff*dcos(theta)+yoff*dsin(theta)) shifty = -(yoff*dcos(theta)-xoff*dsin(theta)) endif do i=1,icount xproj(i) = xproj(i) + shiftx yproj(i) = yproj(i) + shifty enddo call compute_shifted_coords return end subroutine compute_shifted_coords implicit real*8 (a-h,o-z) include 'moe.common.blocks' parameter (raddeg= 57.29577951) real*8 dy_arcsec, delta_y real*8 ra_shift, dec_shift, ra_shift2, dec_shift2 if (.not.coords) return c compute unprojected ra and dec shifts in arc seconds ra_shift = shiftx*dcos(theta) - shifty*dsin(theta) dec_shift = shifty*dcos(theta) + shiftx*dsin(theta) c now compute shift in coordinates due to mask field cut 0.700 inches c in +y direction (so coordinates need to go in -y direction) c delta_y = 0.700 ! old ldss2 mask y zero-point shift delta_y = 0.000 dy_arcsec = delta_y * 25.4 * scale ra_shift2 = -dy_arcsec * dsin(theta) dec_shift2 = dy_arcsec * dcos(theta) c convert to degrees call hms_to_deg(rah,ram,ras,ra_deg) call dms_to_deg(decd,decm,decs,dec_sign,dec_deg) c compute new ra and dec dec_proj = dec_deg - dec_shift/3600. - dec_shift2/3600. ra_proj = ra_deg - ra_shift/(3600.*dcos(dec_deg/raddeg)) % - ra_shift2/(3600.*dcos(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,dec_sign,decdproj,decmproj,decsproj) return end subroutine hms_to_deg(rah,ram,ras,ra_deg) implicit real*8 (a-h,o-z) real*8 rah,ram,ras,ra_deg ra_deg = 15.*(rah + ram/60. + ras/3600.) return end subroutine dms_to_deg(decd,decm,decs,dec_sign,dec_deg) implicit real*8 (a-h,o-z) real*8 decd,decm,decs,dec_sign,dec_deg dec_deg = decd + decm/60. + decs/3600. dec_deg = dec_sign * dec_deg return end subroutine deg_to_hms(ra_deg,rah,ram,ras) implicit real*8 (a-h,o-z) real*8 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,sign,decd,decm,decs) implicit real*8 (a-h,o-z) real*8 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. return end subroutine remove_out_of_field_objects(idone) implicit real*8 (a-h,o-z) include 'moe.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 (priority(kcount(i)).ge.0.0) then c reject object if spectrum does not fall on CCD if (xspecr_ccd(i) .lt. (-xccd_size/2.0) + .or. xspecb_ccd(i) .gt. ( xccd_size/2.0) + .or. yproj(i) .lt. (-0.5*ysize + yedge) + .or. yproj(i) .gt. (0.5* ysize - yedge) ) then ifull = kcount(i) iflag(ifull) =1 idone =0 endif else c reject align star if image not on ccd if ( xmask_ccd(i) - 0.5*(10.0+starwidth/wfpix) .lt. (-xccd_size/2.0) ) then ! demand a 10 pixel x-edge write(6,'(41h ****WARNING: alignment star out of field)') write(6,'(15h object number ,i10,/)')inumobj(kcount(i)) iflag(kcount(i)) =1 idone =0 endif endif ! .... done with field of view check for this object c old code... c if(auto_priority) then c ifull = kcount(i) c if (priority(ifull).lt.0.0 .and. iwarn.ne.1) then c write(6,'(41h ****WARNING: alignment star out of field)') c write(6,'(15h object number ,i10,/)')inumobj(ifull) c endif c iflag(ifull) =1 c idone =0 c endif 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 c real*8 function xmax(y) cc cc compute the maximum x extent from the center (+ve number) of the LDSS mask cc c c include 'moe.common.blocks' c cc xmax = dsqrt( (88.0*scale)**2 - abs(y)**2 ) cc if (xmax.lt.195.) xmax = 195. cc cc the following change to partially account for restricted field cc found by Swara cc c xmax = dsqrt( (80.0*scale)**2 - abs(y)**2 ) c if (xmax.lt.162.) xmax = 162. c c return c end real*8 function xmax(y) c c compute the maximum x extent from the center (+ve number) of the LDSS mask c implicit real*8 (a-h,o-z) include 'moe.common.blocks' c xmax = dsqrt( (88.0*scale)**2 - abs(y)**2 ) c if (xmax.lt.195.) xmax = 195. c c the following change to partially account for restricted field c found by Swara c xmax = dsqrt( (80.0*scale)**2 - abs(y)**2 ) if (xmax.gt.0.5*xsize) xmax = 0.5*xsize return end subroutine compute_slit_ends(idone) implicit real*8 (a-h,o-z) include 'moe.common.blocks' c 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... xstart(i) = xproj(i) - slitmin/2.0 xend(i) = xproj(i) + slitmin/2.0 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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' integer i,ilast,inext ilast = 0 inext = 0 c find closest high side object xmin = xccd_size do j = i+1,icount c if (iflag(kcount(j)).ne.1) then ydist = dabs(yproj(j)-yproj(i)) xdist = dabs(xspecr_ccd(j)-xspecr_ccd(i)) c xdist = dabs(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 = xccd_xsize do j = 1,i-1 if (iflag(kcount(j)).ne.1) then ydist = dabs(yproj(j)-yproj(i)) xdist = dabs(xspecr_ccd(j)-xspecr_ccd(i)) c xdist = dabs(xproj(j)-xproj(i)) if (ydist .le. ystack .and. xdist.le.xmin) then ilast = j xmin = xdist endif endif enddo return end subroutine compute_yends(i) c c compute ystart and yend given the mask PA, and PA of i-th slit c implicit real*8 (a-h,o-z) integer i c *** DOES NOT YET WORK WITH MOE include 'moe.common.blocks' parameter (raddeg= 57.29577951) if (tilted_slits) then alpha = (thinrot-180.-paproj(i) )/raddeg ystart(i) = yproj(i) - (xproj(i) - xstart(i)) * dtan(alpha) yend(i) = yproj(i) + (xend(i) - xproj(i)) * dtan(alpha) else ystart(i) = yproj(i) yend(i) = yproj(i) endif return end subroutine delete_alignment_clashes(idone) implicit real*8 (a-h,o-z) include 'moe.common.blocks' 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 ( xspecr_ccd(i).lt.xspecb_ccd(j) % .and. projpriority(j).ge.0.0 % .and. dabs(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 ( xspecr_ccd(j).lt.xspecb_ccd(i) % .and. projpriority(j).ge.0.0 % .and. dabs(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) implicit real*8 (a-h,o-z) include 'moe.common.blocks' real*8 dx integer idone,ifull,i,inext,ilast 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) then if( xspecr_ccd(inext).lt.xspecb_ccd(i) ) then ! not enough dx--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 add_objects_back implicit real*8 (a-h,o-z) include 'moe.common.blocks' c c A routine to fill available space between accepted slits by trying to add previously c rejected objects back into the mask. This is done in reverse order, from high to low x, c in order to catch objects which can fit, but which clashed with a later deleted slit in c the first pass. c do i = icount,2 if (iflag(kcount(i)).eq.1) then call get_adjacent_objects(i,ilast,inext) if (xspecb_ccd(i).lt.xspecr_ccd(inext) .and. % xspecr_ccd(i).gt.xspecb_ccd(ilast)) then iflag(kcount(i)) = 0 elseif (xspecb_ccd(i).gt.xspecr_ccd(inext)) then if (projpriority(i).lt.projpriority(inext)) then iflag(kcount(i)) = 0 iflag(kcount(inext)) = 1 endif elseif (xspecr_ccd(i).lt.xspecb_ccd(ilast)) then if (projpriority(i).lt.projpriority(ilast)) then iflag(kcount(i)) = 0 iflag(kcount(ilast)) = 1 endif endif endif enddo return end subroutine delete_hole_clashes(idone) implicit real*8 (a-h,o-z) include 'moe.common.blocks' c parameter (wfpix = 0.77402) ! the wfccd pixel size in arcsec real*8 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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' integer idone,i,ifull real*8 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) implicit real*8 (a-h,o-z) 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 implicit real*8 (a-h,o-z) include 'moe.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,sign,line,error) i1 = i2 call get_next_number(i1,i2,ram,sign,line,error) i1 = i2 call get_next_number(i1,i2,ras,sign,line,error) c read declination string from line i1 = i2 call get_next_number(i1,i2,decd,sign,line,error) i1 = i2 call get_next_number(i1,i2,decm,sign,line,error) i1 = i2 call get_next_number(i1,i2,decs,sign,line,error) if (.not.error) then write(6,'(a)')'Enter Epoch' read(5,*) epoch coords = .true. return endif enddo return end subroutine get_ra(line) implicit real*8 (a-h,o-z) include 'moe.common.blocks' logical error character*80 line if (line.eq.' ') return error = .false. c read RA string from line i1 = 0 call get_next_number(i1,i2,rah,sign,line,error) i1 = i2 call get_next_number(i1,i2,ram,sign,line,error) i1 = i2 call get_next_number(i1,i2,ras,sign,line,error) ra_as = nint(rah)*3600.0 + nint(ram)*60.0 + ras if (error) then write(6,'(a)')'Error in RA format' stop endif coords = .true. return end subroutine get_dec(line) implicit real*8 (a-h,o-z) include 'moe.common.blocks' logical error character*80 line if (line.eq.' ') return error = .false. c read declination string from line i1 = 0 call get_next_number(i1,i2,decd,dec_sign,line,error) i1 = i2 call get_next_number(i1,i2,decm,sign,line,error) i1 = i2 call get_next_number(i1,i2,decs,sign,line,error) dec_as = nint(decd)*3600.0 + nint(decm)*60.0 + decs dec_as = dec_as*dec_sign if (error) then write(6,'(a)')'Error in DEC format' stop endif coords = .true. return end subroutine get_next_number(i1,i2,anum,sign,line,error) c c find next string delimted by a space or colon and test that it is a number c implicit real*8 (a-h,o-z) character*80 line logical error,is_a_number if (error) return sign = +1.0 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 anum = abs(anum) i2 = ic else error = .true. endif if (index(line(ilast+1:ic-1),'-').ne.0) sign = -1.0 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) implicit real*8 (a-h,o-z) 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) implicit real*8 (a-h,o-z) 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 implicit real*8 (a-h,o-z) include 'moe.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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' c 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 (.not.cosmic) 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 (.not.cosmic) 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,*)asperunit 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 (.not.cosmic) then write(6,'(52h Enter 1 to modify default field size limits, else 0)') write(6,'(49h default values are 510 arcs along slits [x] and )') write(6,'(41h 330 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 implicit real*8 (a-h,o-z) include 'moe.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_spectrum_file implicit real*8 (a-h,o-z) include 'moe.common.blocks' logical color real*4 x1,y1,x2,y2,xsize4,ysize4 c c plot the current working set of vectors--ioption is initialized to zero c xsize4 = 8342.0 ! xpixel width incl. 3x50 pixel spaces ysize4 = 8242.0 ! ypixel height incl. 50 pixel space color = .false. call initialize_photo_plot(xsize4,ysize4,x1,y1,x2,y2) call make_moe_ccd_outline call make_plot_title(x1,y2) call plot_spectral_footprints(color) call plot_moe_spectrum_numbers(color) call pgend return end subroutine make_postscript_file implicit real*8 (a-h,o-z) include 'moe.common.blocks' logical color real*4 x1,y1,x2,y2,xsize4,ysize4 c c plot the current working set of vectors--ioption is initialized to zero c xsize4 = 8342.0 ! xpixel width including 3x50 pixel spaces ysize4 = 8242.0 ! ypixel height including a 50 pixel space color = .false. call initialize_photo_plot(xsize4,ysize4,x1,y1,x2,y2) call make_moe_ccd_outline call make_plot_title(x1,y2) call plot_slits(color) call plot_numbers(color) call pgend return end subroutine rename_photomask_pgplot_file implicit real*8 (a-h,o-z) character*67 photomask,line*80 write(6,'(a)')'Enter name of postscript file' read(5,'(a67)')photomask line(:13) = 'mv pgplot.ps ' line(14:) = photomask call system(line) return end subroutine make_xwindow_plot implicit real*8 (a-h,o-z) real*8 xsize,ysize real*4 xsize4,ysize4,x1,y1,x2,y2 logical color xsize4 = 8342.0 ! xpixel width including 3x50 pixel spaces ysize4 = 8242.0 ! ypixel height including a 50 pixel space color = .true. call initialize_xwindow_plot(xsize4,ysize4,x1,y1,x2,y2) call make_moe_ccd_outline call make_plot_title(x1,y2) call plot_slits(color) call plot_numbers(color) call pgend return end subroutine make_xwindow_spectrum_plot implicit real*8 (a-h,o-z) real*8 xsize,ysize real*4 xsize4,ysize4,x1,y1,x2,y2 logical color xsize4 = 8342.0 ! xpixel width including 3x50 pixel spaces ysize4 = 8242.0 ! ypixel height including a 50 pixel space color = .true. call initialize_xwindow_plot(xsize4,ysize4,x1,y1,x2,y2) call make_moe_ccd_outline call make_plot_title(x1,y2) call plot_spectral_footprints(color) call plot_moe_spectrum_numbers(color) call pgend return end subroutine initialize_xwindow_plot(xsize4,ysize4,x1,y1,x2,y2) implicit real*8 (a-h,o-z) real*4 xsize4,ysize4,x1,y1,x2,y2 call pgbegin(0,'/xw',1,1) ! open the xwindow x1 = -0.5*xsize4 - 8.5 x2 = +0.5*xsize4 + 8.5 y1 = -0.5*ysize4 - 8.5 y2 = +0.5*ysize4 + 8.5 just = 1 iaxis = -1 call pgslw(1) c x1 = -x1 ! for some reason I need to flip x-direction c x2 = -x2 ! ... maybe not call pgenv(x1,x2,y1,y2,just,iaxis) return end subroutine make_ldss2_mask_outline (ysize4) implicit real*8 (a-h,o-z) real*4 x(1001),y(1001),ysize4,dy dy = ysize4/499.0 y(1) = -0.5*ysize4 x(1) = xmax(dble(0.5*ysize4)) y(1000) = -0.5*ysize4 x(1000) = -xmax(dble(0.5*ysize4)) y(1001) = -0.5*ysize4 x(1001) = xmax(dble(0.5*ysize4)) do i = 2,500 y(i)= y(i-1) + dy y(1000-i+1) = y(i) x(i)= xmax(dble(y(i))) x(1000-i+1) = -xmax(dble(y(i))) enddo call pgline(1001,x,y) return end subroutine make_moe_ccd_outline c this routine draws the outline of the CCDs in pixels real*4 x(5),y(5) c ccd 1? x(1) = -4171.0 y(1) = -4121.0 x(2) = -4171.0 y(2) = -25.0 x(3) = -2123.0 y(3) = -25.0 x(4) = -2123.0 y(4) = -4121.0 x(5) = -4171.0 y(5) = -4121.0 call pgline(5,x,y) c ccd 2? x(1) = -2073.0 y(1) = -4121.0 x(2) = -2073.0 y(2) = -25.0 x(3) = -25.0 y(3) = -25.0 x(4) = -25.0 y(4) = -4121.0 x(5) = -2073.0 y(5) = -4121.0 call pgline(5,x,y) c ccd 3? x(1) = 25.0 y(1) = -4121.0 x(2) = 25.0 y(2) = -25.0 x(3) = 2073.0 y(3) = -25.0 x(4) = 2073.0 y(4) = -4121.0 x(5) = 25.0 y(5) = -4121.0 call pgline(5,x,y) c ccd 4? x(1) = 2123.0 y(1) = -4121.0 x(2) = 2123.0 y(2) = -25.0 x(3) = 4171.0 y(3) = -25.0 x(4) = 4171.0 y(4) = -4121.0 x(5) = 2123.0 y(5) = -4121.0 call pgline(5,x,y) c ccd 5'? x(1) = -4171.0 y(1) = 25.0 x(2) = -4171.0 y(2) = 4121.0 x(3) = -2123.0 y(3) = 4121.0 x(4) = -2123.0 y(4) = 25.0 x(5) = -4171.0 y(5) = 25.0 call pgline(5,x,y) c ccd 6'? x(1) = -2073.0 y(1) = 25.0 x(2) = -2073.0 y(2) = 4121.0 x(3) = -25.0 y(3) = 4121.0 x(4) = -25.0 y(4) = 25.0 x(5) = -2073.0 y(5) = 25.0 call pgline(5,x,y) c ccd 7'? x(1) = 25.0 y(1) = 25.0 x(2) = 25.0 y(2) = 4121.0 x(3) = 2073.0 y(3) = 4121.0 x(4) = 2073.0 y(4) = 25.0 x(5) = 25.0 y(5) = 25.0 call pgline(5,x,y) c ccd 8'? x(1) = 2123.0 y(1) = 25.0 x(2) = 2123.0 y(2) = 4121.0 x(3) = 4171.0 y(3) = 4121.0 x(4) = 4171.0 y(4) = 25.0 x(5) = 2123.0 y(5) = 25.0 call pgline(5,x,y) return end subroutine plot_moe_spectrum_numbers(color) implicit real*8 (a-h,o-z) include 'moe.common.blocks' real*4 x1,y1,x(2),y(2) logical color character*6 text call pgsch(1.0) if (color) then call pgsci(1) call pgscr(1,1.00,1.00,0.00) endif do i=1,icount text=' ' write(text,'(i6)') iobjproj(i) x1 = (xspecr_ccd(i)+xspecb_ccd(i))/2.0 c y1 = (yproj(i)/wfpix - 200.0) y1 = ymask_ccd(i) - 200.0 c switch coords for nod and shuffle orientation if (nod_shuffle) then xdum = x1 ydum = y1 c x1 = -ydum c x1 = ydum c c changed a guess whether it is the spectrum or the number that needs to switch sign c AMcW 17/Nov/2005 c y1 = -xdum ! this was the old y1 = -xdum - 80.0 x1 = -ydum +200.0 c endif c if (iobjproj(i) .ge. 10000) then call pgptext(x1,y1,0.,0.5,text(1:6)) elseif (iobjproj(i) .ge. 1000) then call pgptext(x1,y1,0.,0.5,text(2:6)) else if (iobjproj(i) .ge. 100) then call pgptext(x1,y1,0.,0.5,text(3:6)) else if (iobjproj(i) .ge. 10) then call pgptext(x1,y1,0.,0.5,text(4:6)) else if(iobjproj(i).ge.0) then call pgptext(x1,y1,0.,0.5,text(5:6)) else call pgptext(x1,y1,0.,0.5,text) endif x(1) = xspecb_ccd(i) x(2) = xspecr_ccd(i) y(1) = ymask_ccd(i) y(2) = ymask_ccd(i) if (nod_shuffle) then xdum1 = x(1) xdum2 = x(2) ydum1 = y(1) ydum2 = y(2) x(1) = -ydum1 x(2) = -ydum2 c y(1) = -xdum1 works, but check that it is valid c y(2) = -xdum2 y(1) = -xdum1 y(2) = -xdum2 endif call pgline(2,x,y) enddo call pgsch(1.0) return end subroutine plot_numbers(color) implicit real*8 (a-h,o-z) include 'moe.common.blocks' real*4 x1,y1 logical color character*6 text call pgsch(1.0) c required because pgplot files don't seem to produce color postscript files if (color) then call pgsci(1) call pgscr(1,1.00,1.00,0.00) endif do i=1,icount text=' ' write(text,'(i6)') iobjproj(i) c x1 = xproj(i)/wfpix c y1 = (yproj(i)-30.)/wfpix x1 = xmask_ccd(i) y1 = ymask_ccd(i) - 300.0 y1 = ymask_ccd(i) c switch coords for nod and shuffle orientation if (nod_shuffle) then xdum = x1 ydum = y1 x1 = -ydum - 80. y1 = -xdum + 40 endif c if (iobjproj(i) .ge. 10000) then call pgptext(x1,y1,0.,0.5,text(1:6)) elseif (iobjproj(i) .ge. 1000) then call pgptext(x1,y1,0.,0.5,text(2:6)) else if (iobjproj(i) .ge. 100) then call pgptext(x1,y1,0.,0.5,text(3:6)) else if (iobjproj(i) .ge. 10) then call pgptext(x1,y1,0.,0.5,text(4:6)) else if(iobjproj(i).ge.0) then call pgptext(x1,y1,0.,0.5,text(5:6)) else call pgptext(x1,y1,0.,0.5,text) endif c y1=yproj(i)/wfpix c call pgpoint(1,x1,y1,2) enddo call pgsch(1.0) return end subroutine initialize_photo_plot(xsize4,ysize4,x1,y1,x2,y2) real*4 xsize4,ysize4,x1,y1,x2,y2 call pgbegin(0,'/ps',1,1) ! open the xwindow x1 = -0.5*xsize4 - 100.0 x2 = +0.5*xsize4 + 100.0 y1 = -0.5*ysize4 - 100.0 y2 = +0.5*ysize4 + 100.0 just = 1 iaxis = -1 call pgslw(1) c x1 = -x1 ! for postscript file I don't need to flip origin c x2 = -x2 call pgenv(x1,x2,y1,y2,just,iaxis) return end subroutine make_plot_title(x2,y2) implicit real*8 (a-h,o-z) include 'moe.common.blocks' real*8 rotangle real*4 x2,y2 c character*32 string2 character*50 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+400.0,title) ! -12.0 to +48.0 if(maskpa.gt.359.99999999) maskpa = maskpa - 360.0 c c set estimated rotator offset angle c c 136.2 came rom M.Phillips 29/Oct/2004 c rotangle = maskpa - 136.2 write(string3,155) maskpa, rotangle 155 format('Mask PA = ',f7.2,' Rotator Offset = ',f7.2) c call pgtext(x2-20.,y2+18.0,string3) ! -36.0 to +24.0 call pgtext(x2-20.,y2+220.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 if (dec_sign.lt.0.0) then write(string4,150)irah,iram,dumras,'-', % idecd,idecm,dumdecs,epoch else write(string4,150)irah,iram,dumras,'+', % idecd,idecm,dumdecs,epoch endif else write(string4,152) x0unproj,y0unproj endif 150 format('RA ',i2,1x,i2,1x,f5.2,6x,'DEC ',a1,i2,1x,i2,1x,f4.1,1x, % '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+70.0,string4) ! -48.0 to +12.0 call pgsch(1.0) return end subroutine plot_spectral_footprints(color) c c a routine to plot the outline of the spectra on the IMACS array c include 'moe.common.blocks' real*4 xpts(2),ypts(2) logical color do i = 1,icount c first draw red edge of orders xpts(1) = xspecr_ccd(i) ypts(1) = -4121.0 xpts(2) = xspecr_ccd(i) ypts(2) = 4121.0 if (nod_shuffle) then ypts(1) = -xspecr_ccd(i) xpts(1) = -4151.0 ypts(2) = -xspecr_ccd(i) xpts(2) = 4151.0 endif if (color) then call pgsci(1) call pgscr(1,1.00,0.00,0.00) endif call pgline(2,xpts,ypts) c now draw blue side xpts(1) = xspecb_ccd(i) ypts(1) = 4121.0 xpts(2) = xspecb_ccd(i) ypts(2) = -4121.0 if (nod_shuffle) then ypts(1) = -xspecb_ccd(i) xpts(1) = 4171.0 ypts(2) = -xspecb_ccd(i) xpts(2) = -4171.0 endif if (color) then call pgsci(1) call pgscr(1,0.00,0.00,1.00) endif call pgline(2,xpts,ypts) c return to white color if (color) then call pgsci(1) call pgscr(1,1.00,1.00,1.00) endif enddo return end subroutine plot_slits(color) c now in pixel coordintaes implicit real*8 (a-h,o-z) include 'moe.common.blocks' real*4 xpts(5),ypts(5) real*8 w logical color c required because pgplot files don't seem to produce color postscript files if (color) then call pgsci(1) call pgscr(1,1.00,1.00,0.00) else call pgsci(1) call pgscr(1,0.00,0.00,0.00) endif do i=1,icount c for tilted slits use the projected y-width if (holes .and. projpriority(i).ge.0.0) then c call plot_hole(xproj(i),yproj(i),wproj(i)) !need to change to pixel space else w = wproj(i) w = wslit(kcount(i)) if (ystart(i).ne.yend(i)) then w=w*dsqrt(((ystart(i)-yend(i))/(xstart(i)-xend(i)))**2+1.0) endif xpts(1)=xstart(i)/wfpix c ypts(1)= (ystart(i) - w/2.)/wfpix ypts(1) = ymask_ccd(i) - 0.5*w/wfpix xpts(2)=xstart(i)/wfpix c ypts(2)= (ystart(i) + w/2.)/wfpix ypts(2) = ymask_ccd(i) + 0.5*w/wfpix xpts(3)=xend(i)/wfpix c ypts(3)= (yend(i) + w/2.)/wfpix ypts(3)= ymask_ccd(i) + 0.5*w/wfpix xpts(4)=xend(i)/wfpix c ypts(4)= (yend(i) - w/2.)/wfpix ypts(4)= ymask_ccd(i) - 0.5*w/wfpix xpts(5)=xstart(i)/wfpix c ypts(5)= (ystart(i) - w/2.)/wfpix ypts(5) = ymask_ccd(i) - 0.5*w/wfpix c switch coords for nod and shuffle orientation if (nod_shuffle) then do k = 1,5 xdum = xpts(k) ydum = ypts(k) xpts(k) = -ydum ypts(k) = xdum enddo endif 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 implicit real*8 (a-h,o-z) parameter (raddeg= 57.29577951) real*8 xc,yc,w,r real*4 x(25),y(25) r = w/2. do i = 1,25 alpha = real(i-1)*360./24. x(i) = xc + r*dcos(alpha/raddeg) y(i) = yc + r*dsin(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) implicit real*8 (a-h,o-z) real*8 scale real*4 xgrid(4),ygrid(4) 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) implicit real*8 (a-h,o-z) 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,'(44h 4 to plot spectral footprints on CCD, )') write(6,'(45h 5 type current list to term/outputfile,)') write(6,'(48h 6 to write a postscript slit image file , )') write(6,'(55h 7 to write a postscript spectral footprint file, )') 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,'(40h 11 to plot image of slits on CCD, )') write(6,'(44h 12 to exit program, no further output:,//)') read(5,*)ioption return end subroutine change_slit_pa implicit real*8 (a-h,o-z) include 'moe.common.blocks' parameter (raddeg= 57.29577951) do while (.true.) write(6,'(52h Enter the Slit Position Angle for your observations)') write(6,'(48h acceptable angles are 0.0 < ringangle <= 360.0 )') write(6,'(51h 90 gives N/S slits; 180 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.0.0 .and. thinrot.le.360.0) then thetad = thinrot maskpa = thetad - 90.0 if (maskpa.lt.0.0) maskpa = maskpa + 360.0 theta = thetad/raddeg return endif write(6,'(36h illegal Slit Position Angle entered)') enddo return end subroutine print_current_list c c prints object information to terminal, and to a file if required c implicit real*8 (a-h,o-z) 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 implicit real*8 (a-h,o-z) include 'moe.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 Priority Xspec(red) Xspec(blue) ' write(iunit,'(a)') % ' (arcs) (arcs) (arcs) (arcs) (pix) (pix) ' write(iunit,'(a)') ' ' do j = 1,icount write(iunit,'(i6,10f9.2)') iobjproj(j),xproj(j), % yproj(j),wproj(j), + lproj(j),projpriority(j),xspecr_ccd(j),xspecb_ccd(j), + xmask_ccd(j),ymask_ccd(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(i6,1x,4(f9.2,1x)) 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 implicit real*8 (a-h,o-z) 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) implicit real*8 (a-h,o-z) 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 compute_slit_xymm(xobj,yobj,odec,dec_as,fl,xmm,ymm) implicit real*8 (a-h,o-z) parameter (raddeg= 57.29577951) real*8 fdistort,xobj,yobj,odec,dec_as,fl,xmm,ymm decr = odec/raddeg dec0r = dec_as/(3600.0*raddeg) call compute_xi_eta(xobj,decr,dec0r,fl,xi,eta) call compute_xy_distortions(xobj,yobj,fdistort) c informed guess xmm = -xi * fdistort ymm = -eta * fdistort return end subroutine compute_spectrum_x_positions c c a routine to compute the x-position of the spectra on the ccd plane, with c pixel coordinates which include the 50 pixel gaps between each chip. c implicit real*8 (a-h,o-z) parameter (raddeg= 57.29577951) include 'moe.common.blocks' real*8 fdistort,xpyp,xred,roffset,boffset do i = 1,icount xo = xproj(i)/scale yo = yproj(i)/scale xmask_ccd(i) = -xo * 26.066 ! I need a variable for scale pix/mm ymask_ccd(i) = yo * 26.066 c ymask_ccd(i) = -yo * 26.066 ! -ve for CCD mask image flip ?? c ymask_ccd(i) = 4121.0 + ymask_ccd(i) c xmask_ccd(i) = 4171.0 + xmask_ccd(i) dx = -90.1 + 0.1682*(ymask_ccd(i)+4121.0) ! Tilt effect dx = 102.0 + 0.1662*ymask_ccd(i) ! Tilt effect AMcW 17 Nov 2005 xred = xmask_ccd(i) - dx ! center of red-most slit without blocking filter call xdisp_bounds(xred,roffset,boffset) xspecr_ccd(i) = xred + roffset - 25.0 ! x-disp red bound xspecb_ccd(i) = xred + boffset + 25.0 ! x-disp blue bound c c Below I've added an extra 5 pixels buffer because for the Dec/05 run the parameters have not c yet been empirically determined. I need to update this, and the tilt later. AMcW 17/Nov/2005. c xspecr_ccd(i) = xred + roffset - 30.0 ! x-disp red bound AMcW 17/Nov/05 xspecb_ccd(i) = xred + boffset + 30.0 ! x-disp blue bound enddo return end subroutine compute_xi_eta(xobj,decr,dec0r,fl,xi,eta) implicit real*8 (a-h,o-z) parameter (raddeg= 57.29577951) c dec0r is the declination of the plate center in radians c decr is the object declination in radians c x is the RA difference in radians x = xobj/(raddeg*3600.d0) / dcos(decr) xdum = dsin(dec0r) + (1.0/dtan(decr)) * dcos(dec0r)* dcos(x) xi = fl * dsin(x)*(1.0/dtan(decr)) / xdum eta = fl * ( dcos(dec0r) - (1.0/dtan(decr))*dsin(dec0r)*dcos(x)) / xdum return end subroutine compute_xy_distortions(xobj,yobj,fdistort) c c A routine to compute the distortion factors for the object c implicit real*8 (a-h,o-z) parameter (raddeg= 57.29577951) c convert to radians x = xobj/(raddeg*3600.d0) y = yobj/(raddeg*3600.d0) c compute approximate radius vector length r = dacos( dcos(x)*dcos(y) ) r2 = r*r c telescope distortion dtel = 1.0d0 + r2 * 2.67982953d+02 + r2**2 * 3.68633337d+05 c temperature distortion temp = 5.0 dtemp = 1.0d0 + 1.55108263d-07 -1.94048247d-08 * temp c total distortion fdistort = dtel * dtemp return end c subroutine compute_slit_xymm(xproj,yproj,fl,xmm,ymm) cc cc A routine to compute the x,y position, in mm, of the object, including cc optical distortions cc c c implicit real*8 (a-h,o-z) c parameter (raddeg= 57.29577951) c cc convert to radians c c x = xproj/(raddeg*3600.d0) c y = yproj/(raddeg*3600.d0) c cc compute radius vector length c c r = dacos( dcos(x)*dcos(y) ) c r2 = r*r c cc r2 = x*x + y*y cc r = dsqrt(r2) c cc telescope distortion c c dtel = 1.0d0 + r2 * 2.67982953d+02 + r2**2 * 3.68633337d+05 c cc temperature distortion c c temp = 5.0 c dtemp = 1.0d0 + 1.55108263d-07 -1.94048247d-08 * temp c cc total distortion c c rnew = r * dtel * dtemp c xmm = fl * rnew * x/r c ymm = fl * rnew * y/r c c return c end subroutine xdisp_bounds(xred,roffset,boffset) implicit real*8 (a-h,o-z) include 'moe.common.blocks' call order_shift(red_limit,delta_r1,delta_r2) call order_shift(blue_limit,delta_b1,delta_b2) roffset = delta_r1 + (delta_r2-delta_r1) * (xred-421.5)/(7324.8-421.5) boffset = delta_b1 + (delta_b2-delta_b1) * (xred-421.5)/(7324.8-421.5) return end subroutine order_shift(wavelength,delta1,delta2) c gives the shift in the cross-dispersion direction from a given wavelength c to the red-most unblocked order for two x-positions in xspec_ccd space implicit real*8 (a-h,o-z) if (wavelength.gt.8300.0) then delta1 = 0.0 delta2 = 0.0 elseif (wavelength.gt.7000.0) then delta1 = 50.0 delta2 = 51.0 elseif (wavelength.gt.6100.0) then delta1 = 100.0 delta2 = 104.0 elseif (wavelength.gt.5300.0) then delta1 = 152.0 delta2 = 157.0 elseif (wavelength.gt.4800.0) then delta1 = 209.0 delta2 = 215.0 elseif (wavelength.gt.4300.0) then delta1 = 272.0 delta2 = 280.0 elseif (wavelength.gt.3900.0) then delta1 = 342.0 delta2 = 351.0 elseif (wavelength.gt.3600.0) then delta1 = 419.0 delta2 = 430.0 else delta1 = 503.0 delta2 = 521.0 endif return end subroutine exclude_all_objects(iflag,inumin) implicit real*8 (a-h,o-z) 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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' 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) implicit real*8 (a-h,o-z) parameter (kmax=4000) ! maximum number of objects in a mask list real*8 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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' c character*64 newfile c write(6,'(a)')' write input list of all objects not on current mask' c write(6,'(a)')' but including alignment stars: [1=y, 0=n]' c read(5,*)ilist c if(ilist.eq.1) then if(exclfile.ne.' ') then c write(6,'(a)')' Enter filename for new list' c read(5,'(a)')newfile c 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) implicit real*8 (a-h,o-z) 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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' real*8 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 scalein = 1.0/(25.4*scale) ! inches per arcsec 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 if (ye(i).gt.0.5*ysize*scalein) then dxdy = (xe(i)-xs(i))/(ye(i)-ys(i)) 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 dxdy = (xe(i)-xs(i))/(ye(i)-ys(i)) dys = ys(i) + 0.5*ysize * scalein xs(i) = xs(i) + dys * dxdy ys(i) = - 0.5*ysize * scalein endif enddo c c give tool diameter c if (sw074) then td = 0.010 elseif (sw103) then td = 0.014 elseif (sw125) then td = 0.017 elseif (sw147) then td = 0.020 elseif (sw177) then td = 0.024 elseif (sw199) then td = 0.027 else td=0.020 endif c c re-order slits/holes for efficient milling path c call short_path(xs,ys,xe,ye,sw) call gmcode_2p(xs,ys,xe,ye,sw,td,icount,title) return end subroutine short_path(xs,ys,xe,ye,sw) implicit real*8 (a-h,o-z) include 'moe.common.blocks' real*8 xs(kmax),ys(kmax),xe(kmax),ye(kmax),sw(kmax) real*8 xsn(kmax),ysn(kmax),xen(kmax),yen(kmax),swn(kmax) integer iorder(kmax),i c use Numerical Recipes code for optimal order do i = 1,icount iorder(i) = i enddo call anneal(xs,ys,iorder,icount) c put slits/holes in right order do i = 1,icount xsn(i) = xs(iorder(i)) ysn(i) = ys(iorder(i)) xen(i) = xe(iorder(i)) yen(i) = ye(iorder(i)) swn(i) = sw(iorder(i)) enddo do i = 1,icount xs(i) = xsn(i) ys(i) = ysn(i) xe(i) = xen(i) ye(i) = yen(i) sw(i) = swn(i) enddo 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 implicit real*8 (a-h,o-z) include 'moe.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 call shift_xy ! apply x,y offsets call compute_spectrum_x_positions call sort_slits 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 add_objects_back call check_slit_lengths(idone) endif endif enddo write(6,110)iterate, icount, x0unproj, y0unproj ! is this really needed??? 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 implicit real*8 (a-h,o-z) include 'moe.common.blocks' logical numbers c c identification postscript chart c numbers = .true. iunit=41 call make_postscript_spectrum_file call rename_pgplot_file c call produce_gmcode_file call print_smf_file call print_excluded_objects call print_list(iunit) return end subroutine rename_pgplot_file implicit real*8 (a-h,o-z) include 'moe.common.blocks' character line*80 line(:13) = 'mv pgplot.ps ' line(14:) = psfile call system(line) return end subroutine print_smf_file implicit real*8 (a-h,o-z) include 'moe.common.blocks' call open_smf_file call print_position_to_smf call print_smf_holes_and_slits return end subroutine open_smf_file c open SMF output file implicit real*8 (a-h,o-z) include 'moe.common.blocks' character*(64) smf_file,literal smf_file = literal(gmoutput) do i = 1,68 if (smf_file(i:i).eq.' ') then smf_file(i:) = '.SMF' open(unit=90,file=smf_file(:i+3)) return endif enddo return end subroutine print_position_to_smf c write POSITION keyword to header implicit real*8 (a-h,o-z) include 'moe.common.blocks' character*80 string4 c get the hexadecimal coordinae values 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 if (dec_sign.lt.0.0) then write(string4,150)irah,iram,dumras,'-', % idecd,idecm,dumdecs,epoch,maskpa else write(string4,150)irah,iram,dumras,'+', % idecd,idecm,dumdecs,epoch,maskpa endif 150 format('POSITION 0 ',i2,':',i2,':',f6.3,1x,a1,i2,':',i2,':',f5.2,1x, % f6.1,1x,f5.1,' 0.0 0.0') ! *** need to understand last 3 c write POSITION keyword to header ilength = len(string4) write(90,'(a)')string4(:ilength) return end subroutine print_smf_holes_and_slits implicit real*8 (a-h,o-z) include 'moe.common.blocks' character*90 line do i = 1,icount line = ' ' if (projpriority(i).lt.0.0) then line(:5) = 'HOLE ' else line(:5) = 'SLIT ' endif c write star number write(line(6:13),'(i8)')inumobj(kcount(i)) c write celestial coordinates line(15:26) = ora_str(kcount(i)) line(28:39) = odec_str(kcount(i)) c write slit position in mm c write(line(40:57),'(2f9.3)')xmm(i),ymm(i) write(line(40:57),'(2f9.3)')xmmproj(i),ymmproj(i) c write slit length and width in mm sw = wslit(kcount(i))/scale write(line(59:63),'(f5.3)')sw if (projpriority(i).lt.0.0) then sle = 0.5 * wslit(kcount(i))/scale slw = 0.5 * wslit(kcount(i))/scale else sle = 2.5d0/scale slw = 2.5d0/scale endif if (projpriority(i).lt.0.0) then write(line(64:84),'(3h 1 ,2(1x,f5.3),6h 0.0)')sle,slw write(90,'(a84)')line(:84) else write(line(64:81),'(2(1x,f5.3),6h 0.0)')sle,slw write(90,'(a81)')line(:81) endif enddo return end C************************************************************************ C NOTE*** C C Temporary version to take the 0.007 inch endmill to do hole cutting C SUBROUTINE gmcode_2p(xs,ys,xe,ye,sw,td,NS,title) implicit real*8 (a-h,o-z) 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.012, 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 TD is the tool diameter c c INTEGER LINE,NPOINTS,NSWEEPS,NZIGZAGS,I,IFILE parameter(kmax=4000) REAL*8 XS(kmax),YS(kmax),XE(kmax),YE(kmax),SW(kmax),td CHARACTER title*32 c c default milling parameters c ZMETAL = 0.012 TIP = 0.003 ENDMILLD = 0.020 UNCUT = 0.003 c ZSAFE = 0.100 ZSAFE = 0.010 ZSTART = 0.003 !may want to change zstart to 0.005 ZENDMILL = ZMETAL - UNCUT C ...Open output file c single pass mode for all but 0.010 slits, as the 0.010 milling tool is too short. if (td.gt.0.011) then c write header information to GM code file(s) write(20,'(a,a)') '(', title write(20,'(a,f5.3,a)') '( Single pass mode: ',td,' inch small milling 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)') 'F 1.5 E0.1875' ! note the Z motion is 10% of the X-Y motion write(20,'(a,f5.1)') ' Z ', ZSAFE write(20,'(a)') 'G39 X -3.47 Y - 2.25 I 6.94 J 4.50' ! display ldss mask AMcW 02/2001 c write slit information if (ns.gt.kmax) stop ' TOO MANY SLITS' ms = 2 do ks =1, ns call slit(xs(ks),ys(ks),xe(ks),ye(ks),sw(ks),zsafe, % zstart,zmetal,tip,td,ms,zendmill) enddo c move tool out of the way write(20,'(a)') ' Z 0.80' write(20,'(a)') ' X 2.5 Y -1.6' write(20,'(a)') ' M30' else WRITE(20,'(a,a)') '(', title WRITE(20,'(a,f5.3,a)') '( 2nd pass: ',td,' inch small milling tool' WRITE(30,'(a,a)') '(', title WRITE(30,'(a,f5.3,a)') '( 1st pass: ',ENDMILLD,' ENDMILL' 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)') 'F 1.5 E0.1875' ! note the Z motion is 10% of the X-Y motion WRITE(20,'(a,f5.1)') ' Z ', ZSAFE c WRITE(20,'(a)') 'G39 X -2.3 Y - 1.7 I 4.6 J 3.4' !this displays the full slit AMcW 02/2001 c now trying to accomodate the larger slit mask for LDSS 2 WRITE(20,'(a)') 'G39 X -3.47 Y - 2.25 I 6.94 J 4.50' ! AMcW 02/2001 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.5 E0.375' ! note the Z motion is 10% of the X-Y motion WRITE(30,'(a,f5.1)') ' Z ', ZSAFE c WRITE(30,'(a)') 'G39 X -2.3 Y - 1.7 I 4.6 J 3.4' ! this displays the full slit c ....pattern c now trying to accomodate the larger slit mask for LDSS 2 WRITE(30,'(a)') 'G39 X -3.47 Y - 2.25 I 6.93 J 4.50' ! AMcW 02/2001 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,td,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' endif 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,td, % MS, ZENDMILL) implicit real*8 (a-h,o-z) 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. real*8 td 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 c c Ken Clardy pointed out that one does not want to change inputs in a subroutine c 11/2001 c 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.027) 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 if ( (xe-xs).ne.0.0 .or. (ye-ys).ne.0.0 ) then WRITE(30,101) XC,(XE-XS),YC,(YE-YS) endif c ELSE (sw.gt.0.027) then ( Ooops! AMcW 4/2002 ) ELSE 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: c ZUNDER = 2*(SW - TIP) c c This code should enable milling-tool cutting of the slits/holes AMcW 19/May/00 c code for deciding narrow or wide mode for cutting slits c if (sw.le.0.0273) then knw = 1 else knw = 2 endif zunder = zundernom rad = 0.5 * td c 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. c ZINC = ZBELOW - ZSTART2 !... the (negative) incremental Z motion c c if td > 0.011 we are in single pass mode, so move quickly to just above the metal, c otherwise quick move into the trench c if (td.gt.0.011) then write(20,'(a4,f9.5)') ZC,ZSTART ! move to Z start position rapidly ZINC = ZBELOW - ZSTART !... the (negative) incremental Z motion else WRITE(20,'(a4,f9.5)') ZC,ZSTART2 ! move to Z start position rapidly ZINC = ZBELOW - ZSTART2 !... the (negative) incremental Z motion endif c 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 if (dx.ne.0.0 .or. ddy.ne.0.0) then WRITE(20,101) XC,DX,YC,DDY ! move off start position by DX and DDY endif WRITE(20,102) ZC,ZINC ! cut to final zdepth c IF (XCUT.NE.0.0 .OR. (YE-YS-2.*DDY).NE.0.0) THEN ! if not a hole IF (ABS(XCUT).GE.0.000001 .OR. % ABS(YE-YS-2.*DDY).GE.0.000001) THEN ! if not a hole c 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.0*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************************************************************************ c c Now Numerical Recipes code for Travelling Salesman problem c SUBROUTINE anneal(x,y,iorder,ncity) IMPLICIT REAL*8 (A-H,O-Z) INTEGER ncity,iorder(4000) REAL*8 x(4000),y(4000) CU USES irbit1,metrop,ran3,revcst,revers,trncst,trnspt INTEGER i,i1,i2,idec,idum,iseed,j,k,nlimit,nn,nover,nsucc,n(6), *irbit1 REAL*8 de,path,t,tfactr,ran3,alen,x1,x2,y1,y2 LOGICAL ans alen(x1,x2,y1,y2)=dsqrt((x2-x1)**2+(y2-y1)**2) nover=100*ncity nlimit=10*ncity tfactr=0.9 path=0.0 t=0.5 write(*,*)'Wait for solution to traveling salesman problem...' do 11 i=1,ncity-1 i1=iorder(i) i2=iorder(i+1) path=path+alen(x(i1),x(i2),y(i1),y(i2)) 11 continue i1=iorder(ncity) i2=iorder(1) path=path+alen(x(i1),x(i2),y(i1),y(i2)) idum=-1 iseed=111 do 13 j=1,100 nsucc=0 do 12 k=1,nover 1 n(1)=1+int(ncity*ran3(idum)) n(2)=1+int((ncity-1)*ran3(idum)) if (n(2).ge.n(1)) n(2)=n(2)+1 nn=1+mod((n(1)-n(2)+ncity-1),ncity) if (nn.lt.3) goto 1 idec=irbit1(iseed) if (idec.eq.0) then n(3)=n(2)+int(abs(nn-2)*ran3(idum))+1 n(3)=1+mod(n(3)-1,ncity) call trncst(x,y,iorder,ncity,n,de) call metrop(de,t,ans) if (ans) then nsucc=nsucc+1 path=path+de call trnspt(iorder,ncity,n) endif else call revcst(x,y,iorder,ncity,n,de) call metrop(de,t,ans) if (ans) then nsucc=nsucc+1 path=path+de call revers(iorder,ncity,n) endif endif if (nsucc.ge.nlimit) goto 2 12 continue 2 continue c 2 write(*,*) c write(*,*) 'T =',t,' Path Length =',path c write(*,*) 'Successful Moves: ',nsucc t=t*tfactr if (nsucc.eq.0) return 13 continue return END SUBROUTINE revcst(x,y,iorder,ncity,n,de) IMPLICIT REAL*8 (A-H,O-Z) INTEGER ncity,iorder(4000),n(6) REAL*8 de,x(4000),y(4000) INTEGER ii,j REAL*8 alen,xx(4),yy(4),x1,x2,y1,y2 alen(x1,x2,y1,y2)=dsqrt((x2-x1)**2+(y2-y1)**2) n(3)=1+mod((n(1)+ncity-2),ncity) n(4)=1+mod(n(2),ncity) do 11 j=1,4 ii=iorder(n(j)) xx(j)=x(ii) yy(j)=y(ii) 11 continue de=-alen(xx(1),xx(3),yy(1),yy(3))-alen(xx(2),xx(4),yy(2),yy(4))+ *alen(xx(1),xx(4),yy(1),yy(4))+alen(xx(2),xx(3),yy(2),yy(3)) return END SUBROUTINE revers(iorder,ncity,n) IMPLICIT REAL*8 (A-H,O-Z) INTEGER ncity,iorder(4000),n(6) INTEGER itmp,j,k,l,nn nn=(1+mod(n(2)-n(1)+ncity,ncity))/2 do 11 j=1,nn k=1+mod((n(1)+j-2),ncity) l=1+mod((n(2)-j+ncity),ncity) itmp=iorder(k) iorder(k)=iorder(l) iorder(l)=itmp 11 continue return END SUBROUTINE trncst(x,y,iorder, ncity,n,de) IMPLICIT REAL*8 (A-H,O-Z) INTEGER ncity,iorder(4000),n(6) REAL*8 de,x(4000),y(4000) INTEGER ii,j REAL*8 xx(6),yy(6),alen,x1,x2,y1,y2 alen(x1,x2,y1,y2)=dsqrt((x2-x1)**2+(y2-y1)**2) n(4)=1+mod(n(3),ncity) n(5)=1+mod((n(1)+ncity-2),ncity) n(6)=1+mod(n(2),ncity) do 11 j=1,6 ii=iorder(n(j)) xx(j)=x(ii) yy(j)=y(ii) 11 continue de=-alen(xx(2),xx(6),yy(2),yy(6))-alen(xx(1),xx(5),yy(1), *yy(5))-alen(xx(3),xx(4),yy(3),yy(4))+alen(xx(1),xx(3),yy(1), *yy(3))+alen(xx(2),xx(4),yy(2),yy(4))+alen(xx(5),xx(6),yy(5),yy(6)) return END SUBROUTINE trnspt(iorder,ncity,n) IMPLICIT REAL*8 (A-H,O-Z) INTEGER ncity,iorder(4000),n(6),MXCITY PARAMETER (MXCITY=4000) INTEGER j,jj,m1,m2,m3,nn,jorder(MXCITY) m1=1+mod((n(2)-n(1)+ncity),ncity) m2=1+mod((n(5)-n(4)+ncity),ncity) m3=1+mod((n(3)-n(6)+ncity),ncity) nn=1 do 11 j=1,m1 jj=1+mod((j+n(1)-2),ncity) jorder(nn)=iorder(jj) nn=nn+1 11 continue do 12 j=1,m2 jj=1+mod((j+n(4)-2),ncity) jorder(nn)=iorder(jj) nn=nn+1 12 continue do 13 j=1,m3 jj=1+mod((j+n(6)-2),ncity) jorder(nn)=iorder(jj) nn=nn+1 13 continue do 14 j=1,ncity iorder(j)=jorder(j) 14 continue return END SUBROUTINE metrop(de,t,ans) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 de,t LOGICAL ans CU USES ran3 INTEGER jdum REAL*8 ran3 SAVE jdum DATA jdum /1/ ans=(de.lt.0.0).or.(ran3(jdum).lt.exp(-de/t)) return END FUNCTION ran3(idum) IMPLICIT REAL*8 (A-H,O-Z) INTEGER idum INTEGER MBIG,MSEED,MZ C REAL*8 MBIG,MSEED,MZ REAL*8 ran3,FAC PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG) C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG) INTEGER i,iff,ii,inext,inextp,k INTEGER mj,mk,ma(55) C REAL*8 mj,mk,ma(55) SAVE iff,inext,inextp,ma DATA iff /0/ if(idum.lt.0.or.iff.eq.0)then iff=1 mj=MSEED-iabs(idum) mj=mod(mj,MBIG) ma(55)=mj mk=1 do 11 i=1,54 ii=mod(21*i,55) ma(ii)=mk mk=mj-mk if(mk.lt.MZ)mk=mk+MBIG mj=ma(ii) 11 continue do 13 k=1,4 do 12 i=1,55 ma(i)=ma(i)-ma(1+mod(i+30,55)) if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG 12 continue 13 continue inext=0 inextp=31 idum=1 endif inext=inext+1 if(inext.eq.56)inext=1 inextp=inextp+1 if(inextp.eq.56)inextp=1 mj=ma(inext)-ma(inextp) if(mj.lt.MZ)mj=mj+MBIG ma(inext)=mj ran3=mj*FAC return END FUNCTION irbit1(iseed) IMPLICIT REAL*8 (A-H,O-Z) INTEGER irbit1,iseed,IB1,IB2,IB5,IB18 PARAMETER (IB1=1,IB2=2,IB5=16,IB18=131072) LOGICAL newbit newbit=iand(iseed,IB18).ne.0 if(iand(iseed,IB5).ne.0)newbit=.not.newbit if(iand(iseed,IB2).ne.0)newbit=.not.newbit if(iand(iseed,IB1).ne.0)newbit=.not.newbit irbit1=0 iseed=iand(ishft(iseed,1),not(IB1)) if(newbit)then irbit1=1 iseed=ior(iseed,IB1) endif return END