program map CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C This program will read either zcat.2000.dat or gals2000.dat C C and search the catalogue for entries matching the C C desired input parameters, the program will plot those C C objects on the sky using either a linear or an Aitoff C C projection. Different point styles or different colors C C can be used to signify different velocity intervals. C C C C To compile the program and link to the SuperMongo C C libraries: C C C C f77 map.f -l X -L /opt/lib -l plotsub -l devices -l utils C C C C C C Written by Jeff Mader C C C C Last updated 06/1/02 to zcat.2000 format JPH C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Aitoff Projection C C C C RAin=RAin-RAmid # if not centered around 12h # C C RAin=RAin*15.0 # change to decimal degrees # C C costhe=cosd(DECin) C C w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RAin/2))) C C RAout=2*w*costhe*sind(RAin/2) C C RAout=RAout/15.0 # change back to decimal hours # C C RAout=RAout+RAmid # if not centered around 12h # C C DECout=w*sind(DECin) C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC implicit real*8(a-h,o-z) integer i, j, k, label integer rah, ram, decd, decm, decs integer vels,velup, veldown, v(2000000), velin(2000000) integer legendv, colors integer count, num real temp1, minscale, lweight, epoch, dy, decs2 real RA, DEC, tempmaxra, tempmaxdec, tempmindec, tempmiddec real RAmin, RAmax, RAmid, dRA, dDEC, ras, rain real DECmin, DECmax, decin, sizera, sizedec real magup, magdown, mags, magin(2000000) real point1(2000000) real xin(2000000), xin2(2000000), yin(2000000) real yin2(2000000) real x, x1, x2, y1, y2, xmin, xmax, ymin, ymax, dx real ex1(10), expand(2000000), lmin, lmax, lmid double precision PI, costhe, w character*1 sign,magyn,velyn,ans2,ans3,proj,legend,otation,type2 character*2 label1, typein(2000000), type character*3 label2 character*6 label3 character*10 c1(10), color1(2000000) character*17 name character*21 name2 character*30 title, filein character*40 filename CCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Retrieve input from user C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCC print *,' ' print *,' ' type 100 100 format(' ***** ZCAT2000/GALS2000 MAPPING PROGRAM *****') print *,' ' type 500 500 format(' This program will create a map of the sky, in either') type 501 501 format(' linear or aitoff projections, using data from either') type 502 502 format(' zcat2000 or gals2000. The map will contain those zcat') type 503 503 format(' or gals2000 entries within the specified RA, DEC, velocity -') type 504 504 format(' and magnitude ranges. The plotted points will be', - ' scaled') type 505 505 format(' by magnitude, with the minimum scale size being input,' -) type 506 506 format(' and different point styles will represent different') type 507 507 format(' morphological types. Point colors can be used to repre -sent') type 508 508 format(' different velocity ranges.') print *,' ' print *,' ' C type 101 C 101 format(' Input file (zcat format, "z" for zcat.dat or "g" for C - tmassv3.dat): '$) 509 type 101 101 format(' "z" for zcat.dat or "g" for tmassv3.dat, - or "o" for other: '$) accept *,filein print *,' ' if(filein.ne.'z'.and.filein.ne.'g'.and.filein.ne.'o') goto 509 if(filein.eq.'o') type 1110 1110 format(' input other filename: '$) accept *,filename type 102 102 format(' Linear or Aitoff projection (l/a): '$) accept *,proj print *,' ' type 103 103 format(' Enter ZCAT/GALS2000 search parameters-----') print *,' ' type 104 104 format(' RAmin and RAmax (decimal hours): '$) accept *,RAmin,RAmax type 105 105 format(' DECmin and DECmax (decimal degrees): '$) accept *,DECmin,DECmax type 126 126 format(' Input equinox (will also be output equinox): '$) accept *,epoch type 106 106 format(' Enter lower and upper velocity limits: '$) accept *,veldown,velup if(veldown.le.0.or.velup.le.0) then type 107 107 format(' Include those entries with no velocity (y/n): ' -$) accept *,velyn endif type 108 108 format(' Enter lower and upper magnitude limits: '$) accept *,magdown, magup if(magdown.le.0) then type 109 109 format(' Include those entries with no magnitude (y/n): -'$) accept *,magyn endif CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Change RA and DEC into decimal hours and decimal degrees C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC if(RAmax.gt.RAmin) RA=(RAmax+RAmin)/2.0 if(RAmax.lt.RAmin) RA=((RAmin-24)+RAmax)/2.0 DEC=(DECmax+DECmin)/2.0 if(RAmax.gt.RAmin) sizera=(RAmax-RAmin)*15.0 if(RAmax.lt.RAmin) sizera=-((RAmin-24)-RAmax)*15.0 sizedec=DECmax-DECmin CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Read in data from zcat.dat C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC print *,' ' if(filein.eq.'g') type 110 if(filein.eq.'z') type 210 110 format(' Searching xsc.zcat for entries') 210 format(' Searching zcat.2000.dat for entries') print *,' ' i=1 goto 200 996 write(6,*) ' Error opening input file ' close(unit=1) go to 999 200 if(filein.eq.'z') - open(unit=1,file='/home/huchra/z/zcat.2000.dat',status='old', - err=996) if(filein.eq.'g') - open(unit=1,file='/home/jmader/Work/2MASS/gals2000/gals2000.dat', - status='old',err=996) if(filein.eq.'o') - open(unit=1,file=filename,status='old',err=996) num=1 if(filein.eq.'g') then do 50 j=1,7 read(1,300) temp 300 format(A100) 50 continue num=8 endif do 1 j=num,2000000 if(filein.ne.'g') then read(1,202,END=998,ERR=997) name,rah,ram,ras,sign, - decd,decm,decs,mags,velss,type 202 format(A17,2I2,f4.1,A1,3I2,f5.2,f7.0,8X,A2) endif ivtest=velss vtest=float(velss-ivtest) if(vtest.ne.0.0) velss=velss*299725. vels=int(velss) if(filein.eq.'g') then read(1,301,END=998,ERR=997) name2,rah,ram,ras,sign, - decd,decm,decs2,mags,type2,vels,type 301 format(A21,2X,I2,1X,I2,1X,f5.2,1X,A1,I2,1X,I2,1X,f4.1,27X,f6.3,1 -3X,A1,I5,12X,A2) endif if(filein.eq.'g'.and.(type2.eq.'A'.or.type2.eq.'S'.or.type2.eq.'P -')) goto 1 if(mags.lt.magdown.or.mags.gt.magup) goto 1 if(mags.eq.0.and.magyn.eq.'n') goto 1 if(vels.lt.veldown.or.vels.gt.velup) goto 1 if(vels.eq.0.and.velyn.eq.'n') goto 1 RAin=rah+ram/60.0+ras/3600.0 DECin=decd+decm/60.0+decs/3600.0 if(sign.eq.'-') DECin=-DECin C Precess, if necessary if(epoch.ne.2000) then dy=epoch-2000.0 dRA=46.092+(0.0002797*dy)+((20.051-0.0000834*dy)* -sind(RAin*15.0)*tand(DECin)) dRA=dRA*dy/15.0/3600.0 RAin=RAin+dRA if(RAin.lt.0) RAin=RAin+24 if(RAin.ge.24) RAin=RAin-24 dDEC=(20.051-0.0000834*dy)*cosd(RAin*15.0) dDEC=dDEC*dy/3600.0 DECin=DECin+dDEC endif if(RAmin.lt.RAmax) goto 113 if(RAmin.gt.RAmax) goto 115 C C If RAmin < RAmax C 113 if(RAmax.gt.RAmin) then if(RAin.ge.RAmin.and.RAin.le.RAmax) goto 114 goto 1 114 if(DECin.ge.DECmin.and.DECin.le.DECmax) goto 117 goto 1 endif C C If RAmin > RAmax C 115 if(RAmax.lt.RAmin) then if(RAin.le.RAmax.or.RAin.ge.RAmin) goto 116 goto 1 116 if(DECin.ge.DECmin.and.DECin.le.DECmax) goto 118 goto 1 endif 117 xin(i)=RAin yin(i)=DECin magin(i)=mags velin(i)=vels typein(i)=type i=i+1 c testout = type accepted galaxies c type 202, name,rah,ram,ras,sign, c - decd,decm,decs,mags,velss,type goto 1 118 if(RA.le.RAmin.and.RA.le.RAmax) then if(RAin.ge.RAmin) xin(i)=RAin-24 if(RAin.le.RAmax) xin(i)=RAin yin(i)=DECin magin(i)=mags velin(i)=vels typein(i)=type i=i+1 c testout = type accepted galaxies c type 202, name,rah,ram,ras,sign, c - decd,decm,decs,mags,velss,type goto 1 endif if(RA.ge.RAmin.and.RA.ge.RAmax) then if(RAin.ge.RAmin) xin(i)=RAin if(RAin.le.RAmax) xin(i)=RAin+24 yin(i)=DECin magin(i)=mags velin(i)=vels typein(i)=type i=i+1 c testout = type accepted galaxies c type 202, name,rah,ram,ras,sign, c - decd,decm,decs,mags,velss,type endif 1 continue print *,' ' 997 type 119, j 119 format(/' Error in read at ',i5) 998 num=j-1 if(filein.eq.'g') write (*,120) j-8 if(filein.eq.'o') write (*,129) j-1 if(filein.eq.'z') write (*,220) j-1 120 format(' ',i6,' entries in gals2000.dat') 220 format(' ',i6,' entries in zcat.dat') 129 format(' ',i6,' entries in input file') if(filein.eq.'g') write (*,121) i-1 if(filein.eq.'z') write (*,221) i-1 if(filein.eq.'o') write (*,222) i-1 121 format(' 'i6' entries in gals2000.dat matching input criteria -') 221 format(' 'i6' entries in zcat.dat matching input criteria') print *,' ' 222 format(' 'i6' entries in input file matching input criteria') close(unit=1) CCCCCCCCCCCCCCCCCCCCCC C C C Set up Mongo plots C C C CCCCCCCCCCCCCCCCCCCCCC type 122 122 format(' SuperMongo plot setup-----') print *,' ' type 123 123 format(' Plotting Device (1=x11,2=psp5,3=psp3,4=psc1,5=post -script): '$) accept *,iterm if(iterm.ne.1) then type 124 124 format(' Portrait or landscape (p,l): '$) accept *,otation endif type 906 906 format(' What line weight: '$) accept *,lweight type 125 125 format(' Use different colors for velocity intervals (y/n): - '$) accept *,ans2 type 127 127 format(' Use different sizes for magnitude intervals (y/n): - '$) accept *,ans3 if(ans2.eq.'y') then print *," " type 136 136 format(' **** default(b/w), red, green, blue, cyan, magent -a, yellow ****') print *," " type 137 137 format(' How many different colors (7 max): '$) accept *,colors v(1)=veldown ex1(1)=1.0 c1(1)='default' do 5 j=1,colors write (*,138) j 138 format(' Color#',i1,': '$) accept *,c1(j+1) if(j.eq.colors) then write (*,139) velup 139 format(' Velocity limit is ',i6) v(j+1)=velup goto 142 endif type 140 140 format(' Upper velocity limit: '$) accept *,v(j+1) 5 continue endif 142 print *," " type 905 905 format(' Minimum scale size: '$) accept *,minscale type 907 907 format(' Mag for Minimum scale size: '$) accept *,minmag count=0 do 4 k=1,i-1 point1(k)=30.0 if(typein(k).eq.'-7') point1(k)=203.0 if(typein(k).eq.'-6') point1(k)=203.0 if(typein(k).eq.'-5') point1(k)=203.0 if(typein(k).eq.'-4') point1(k)=203.0 if(typein(k).eq.'-3') point1(k)=203.0 if(typein(k).eq.'-2') point1(k)=200.0 if(typein(k).eq.'-1') point1(k)=200.0 if(typein(k).eq.' 0') point1(k)=200.0 if(typein(k).eq.' 1') point1(k)=81.0 if(typein(k).eq.' 2') point1(k)=81.0 if(typein(k).eq.' 3') point1(k)=81.0 if(typein(k).eq.' 4') point1(k)=81.0 if(typein(k).eq.' 5') point1(k)=81.0 if(typein(k).eq.' 6') point1(k)=41.0 if(typein(k).eq.' 7') point1(k)=41.0 if(typein(k).eq.' 8') point1(k)=41.0 if(typein(k).eq.' 9') point1(k)=41.0 if(typein(k).eq.'10') point1(k)=41.0 c pointsizes expand(k)=minscale if(ans3.eq.'n') go to 888 if(magin(k).ge.minmag) go to 888 expand(k) = minscale + 0.25*(minmag-magin(k)) 888 continue if(ans2.eq.'y') then do 6 j=1,colors if(j.eq.1) then if(velin(k).ge.v(j).and.velin(k).le.v(j+1)) color1(k)=c1(j+1) endif if(j.ne.1) then if(velin(k).gt.v(j).and.velin(k).le.v(j+1)) color1(k)=c1(j+1) endif 6 continue endif count=count+1 4 continue if(ans2.eq.'n') goto 146 print *,' ' 146 type 147 147 format(' Title of plot (one word,A30, x for none): '$) accept *,title print *,' ' type 148 148 format(' Would you like a legend added to the plot (y/n): ' -$) accept *,legend if(RAmax.lt.RAmin) RAmin=RAmin-24 temp1=RAmin RAmin=RAmax RAmax=temp1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Call SuperMongo C C plotting area is 170mm x 170mm for portrait C C plotting area is 227mm x 170mm for landscape C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCC C C C Borders C C C CCCCCCCCCCC if(RAmin.lt.RAmax) dx=0.01 if(RAmin.gt.RAmax) dx=-0.01 PI=3.141592653589793238462643 dRA=(RAmax+RAmin) RAmid=dRA/2.0 if(proj.eq.'a') then j=1 DEC=DECmin do 13 x=RAmin,RAmax,dx RA=x-RAmid RA=RA*15.0 costhe=cosd(DEC) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) xin2(j)=2*w*costhe*sind(RA/2) xin2(j)=xin2(j)/15.0 xin2(j)=xin2(j)+RAmid yin2(j)=w*sind(DEC) j=j+1 13 continue RA=RAmax RA=RA-RAmid RA=RA*15 do 14 DEC=DECmin,DECmax,0.01 costhe=cosd(DEC) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) xin2(j)=2*w*costhe*sind(RA/2) xin2(j)=xin2(j)/15.0 xin2(j)=xin2(j)+RAmid yin2(j)=w*sind(DEC) j=j+1 14 continue DEC=DECmax do 15 x=RAmax,RAmin,-dx RA=x-RAmid RA=RA*15 costhe=cosd(DEC) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) xin2(j)=2*w*costhe*sind(RA/2) xin2(j)=xin2(j)/15.0 xin2(j)=xin2(j)+RAmid yin2(j)=w*sind(DEC) j=j+1 15 continue RA=RAmax RA=RAmin RA=RA-RAmid RA=RA*15 do 16 DEC=DECmax,DECmin,-0.01 costhe=cosd(DEC) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) xin2(j)=2*w*costhe*sind(RA/2) xin2(j)=xin2(j)/15.0 xin2(j)=xin2(j)+RAmid yin2(j)=w*sind(DEC) j=j+1 16 continue xmin=xin2(1) xmax=xin2(1) ymin=yin2(1) ymax=yin2(1) tempmiddec=yin2(1) do 99 k=1,j-1,1 if(xin2(k).lt.xmin) xmin=xin2(k) if(xin2(k).gt.xmax) xmax=xin2(k) if(yin2(k).lt.ymin) ymin=yin2(k) if(yin2(k).gt.ymax) ymax=yin2(k) if(abs(yin2(k)-0).lt.(abs(tempmiddec-0))) tempmiddec=yin2(k) 99 continue xmin=xmin-(1.0/170.0) xmax=xmax+(1.0/170.0) ymin=ymin-(1.0/170.0) ymax=ymax+(1.0/170.0) lmin=ymin lmax=ymax endif if(proj.eq.'l') then xmin=RAmax-(1.0/170.0) xmax=RAmin+(1.0/170.0) ymin=DECmin-(1.0/170.0) ymax=DECmax+(1.0/170.0) endif tempmaxra=xmin tempmindec=ymin tempmaxdec=ymax if(xmax.lt.xmin) xmin=xmin-24 if(xmax.gt.xmin) sizera=(xmax-xmin)*15.0 if(xmax.lt.xmin) sizera=-((xmax+24)-xmin)*15.0 sizedec=ymax-ymin if(sizera.eq.sizedec) then if(otation.eq.'l') then xmin=xmin-sizera*(28.5/227.0)/15.0 xmax=xmax+sizera*(28.5/227.0)/15.0 endif if(legend.eq.'y'.and.otation.ne.'l') then xmin=xmin-(18.0*sizera/15.0)/170.0 ymin=ymin-(9.0*sizedec)/170.0 ymax=ymax+(9.0*sizedec)/170.0 endif if(legend.eq.'y'.and.otation.eq.'l') then xmin=xmin-(18.0*sizera/15.0)/227.0 ymin=ymin-(9.0*sizedec)/170.0 ymax=ymax+(9.0*sizedec)/170.0 endif endif if(sizedec.gt.sizera) then if(otation.eq.'p') then xmin=xmin-(1-sizera/sizedec)*(sizera/15.0) xmax=xmax+(1-sizera/sizedec)*(sizera/15.0) endif if(otation.eq.'l') then if((sizera/227.0).lt.(sizedec/170.0)) then xmin=xmin-(sizera/227.0/15.0) -*((227-(sizera*170.0)/sizedec)/2.0) xmax=xmax+(sizera/227.0/15.0) -*((227-(sizera*170.0)/sizedec)/2.0) endif if((sizera/227.0).gt.(sizedec/170.0)) then ymin=ymin-(sizedec/170.0)*((170-(sizedec*227.0)/sizera)/2.0) ymax=ymax+(sizedec/170.0)*((170-(sizedec*227.0)/sizera)/2.0) endif endif if(legend.eq.'y'.and.otation.ne.'l') then xmin=xmin-(18.0*sizera/15.0)/170.0 ymin=ymin-(9.0*sizedec)/170.0 ymax=ymax+(9.0*sizedec)/170.0 endif if(legend.eq.'y'.and.otation.eq.'l') then xmin=xmin-(18.0*sizera/15.0)/227.0 ymin=ymin-(9.0*sizedec)/170.0 ymax=ymax+(9.0*sizedec)/170.0 endif endif if(sizera.gt.sizedec) then if(otation.eq.'p') then ymin=ymin-(1-sizedec/sizera)*(sizedec) ymax=ymax+(1-sizedec/sizera)*(sizedec) endif if(otation.eq.'l') then if((sizera/227.0).lt.(sizedec/170.0)) then xmin=xmin-(sizera/227.0/15.0) -*((227-(sizera*170.0)/sizedec)/2.0) xmax=xmax+(sizera/227.0/15.0) -*((227-(sizera*170.0)/sizedec)/2.0) endif if((sizera/227.0).gt.(sizedec/170.0)) then ymin=ymin-(sizedec/170.0)*((170-(sizedec*227.0)/sizera)/2.0) ymax=ymax+(sizedec/170.0)*((170-(sizedec*227.0)/sizera)/2.0) endif endif if(legend.eq.'y'.and.otation.ne.'l') then xmin=xmin-(18.0*sizera/15.0)/170.0 ymin=ymin-(9.0*sizedec)/170.0 ymax=ymax+(9.0*sizedec)/170.0 endif if(legend.eq.'y'.and.otation.eq.'l') then xmin=xmin-(18.0*sizera/15.0)/227.0 ymin=ymin-(9.0*sizedec)/170.0 ymax=ymax+(9.0*sizedec)/170.0 endif endif call sm_graphics call sm_defvar("TeX_strings","1") if(iterm.eq.1) call sm_device('x11') if(iterm.eq.2.and.otation.eq.'p') call sm_device('psp5') if(iterm.eq.3.and.otation.eq.'p') call sm_device('psp3') if(iterm.eq.4.and.otation.eq.'p') call sm_device('psc1') if(iterm.eq.2.and.otation.eq.'l') call sm_device('psp5la') if(iterm.eq.3.and.otation.eq.'l') call sm_device('psp3la') if(iterm.eq.4.and.otation.eq.'l') call sm_device('psc1la') if(iterm.eq.5) then if(ans2.eq.'n'.and.otation.eq.'p') - call sm_device('postfile map.ps') if(ans2.eq.'n'.and.otation.eq.'l') - call sm_device('postlandfile map.ps') if(ans2.eq.'y'.and.otation.eq.'p') - call sm_device('postfile_colour map.ps') if(ans2.eq.'y'.and.otation.eq.'l') - call sm_device('postlandfile_colour map.ps') endif call sm_limits(xmax,xmin,ymin,ymax) call sm_lweight(lweight) call sm_ctype('default') if(proj.eq.'l') then call sm_relocate(RAmin,DECmin) call sm_draw(RAmin,DECmax) call sm_draw(RAmax,DECmax) call sm_draw(RAmax,DECmin) call sm_draw(RAmin,DECmin) endif if(proj.eq.'a') call sm_conn(xin2,yin2,j-1) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Plot Points (l=linear, a=aitoff) C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC if(ans2.eq.'n') then do 8 j=1, i-1 if(proj.eq.'a') then RA=xin(j)-RAmid RA=RA*15.0 DEC=yin(j) costhe=cosd(DEC) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) xin(j)=2*w*costhe*sind(RA/2) xin(j)=xin(j)/15.0 xin(j)=xin(j)+RAmid yin(j)=w*sind(DEC) endif call sm_ptype(point1(j),1) call sm_expand(expand(j),1) call sm_points(xin(j),yin(j),1) 8 continue endif if(ans2.eq.'y') then do 11 j=1, i-1 if(proj.eq.'a') then RA=xin(j)-RAmid RA=RA*15.0 DEC=yin(j) costhe=cosd(DEC) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) xin(j)=2*w*costhe*sind(RA/2) xin(j)=xin(j)/15.0 xin(j)=xin(j)+RAmid yin(j)=w*sind(DEC) endif call sm_ptype(point1(j),1) call sm_expand(expand(j),1) if(color1(j).eq.'default') call sm_ctype('black') if(color1(j).eq.'red') call sm_ctype('red') if(color1(j).eq.'green') call sm_ctype('green') if(color1(j).eq.'blue') call sm_ctype('blue') if(color1(j).eq.'cyan') call sm_ctype('cyan') if(color1(j).eq.'magenta') call sm_ctype('magenta') if(color1(j).eq.'yellow') call sm_ctype('yellow') if(color1(j).eq.'black') call sm_ctype('black') call sm_points(xin(j),yin(j),1) 11 continue endif call sm_ctype('default') call sm_expand(1.0) CCCCCCCCC C C C Title C C C CCCCCCCCC RA=(RAmax+RAmin)/2.0 DEC=(DECmax+DECmin)/2.0 if(title.ne.'x') then if(proj.eq.'a') then RA=RA-RAmid RA=RA*15.0 costhe=cosd(DECmax) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x1=2*w*costhe*sind(RA/2) x1=x1/15.0 x1=x1+RAmid y1=w*sind(DECmax) endif if(proj.eq.'l') then x1=RA y1=DECmax endif call sm_relocate (x1,(y1+(10.0*sizedec)/170.0)) call sm_putlabel(5,title) endif CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C Create tick marks and label for each axis C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC sizera=RAmax-RAmin dx=DECmax-DECmin sizedec=ymax-ymin if(abs(dx).ge.160) goto 149 if(abs(sizera).le.0.05) dra=0.01 if(abs(sizera).gt.0.05.and.abs(sizera).le.0.75) dra=0.05 if(abs(sizera).gt.0.75.and.abs(sizera).le.1) dra=0.10 if(abs(sizera).gt.1.and.abs(sizera).le.2) dra=0.25 if(abs(sizera).gt.2.and.abs(sizera).le.5) dra=0.50 if(abs(sizera).gt.5.and.abs(sizera).le.10) dra=1.00 if(abs(sizera).gt.10) dra=2.00 if(abs(sizera).eq.24.and.dx.gt.100) dra=6.00 if(abs(sizera).eq.24.and.dx.gt.140) dra=12.00 do 17 x=-24,24,dra if((x.le.RAmax-0.001.and.x.ge.RAmin+0.001).or - .(x.le.RAmin+0.001.and.x.ge.RAmax-0.001)) - then x1=x x2=x y1=DECmax y2=DECmin if(proj.eq.'a') then RA=x-RAmid RA=RA*15.0 costhe=cosd(y1) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x1=2*w*costhe*sind(RA/2) x1=x1/15.0 x1=x1+RAmid y1=w*sind(y1) costhe=cosd(y2) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x2=2*w*costhe*sind(RA/2) x2=x2/15.0 x2=x2+RAmid y2=w*sind(y2) endif call sm_relocate(x1,y1) dx=0.2 if(proj.eq.'a') then do 18 j=1,10 y1=DECmax-((j*dx)*sizedec)/170.0 if(proj.eq.'a') then costhe=cosd(y1) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x1=2*w*costhe*sind(RA/2) x1=x1/15.0 x1=x1+RAmid y1=w*sind(y1) endif call sm_draw(x1,y1) 18 continue endif if(proj.eq.'l') then call sm_draw(x1,y1-2.0*sizedec/170.0) endif call sm_relocate(x2,y2) if(proj.eq.'a') then do 19 j=1,10 y2=DECmin+((j*dx)*sizedec)/170.0 if(proj.eq.'a') then costhe=cosd(y2) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x2=2*w*costhe*sind(RA/2) x2=x2/15.0 x2=x2+RAmid y2=w*sind(y2) endif call sm_draw(x2,y2) 19 continue endif if(proj.eq.'l') then call sm_draw(x2,y2+2.0*sizedec/170.0) endif y2=DECmin-(2.5*sizedec)/170.0 if(proj.eq.'a') then costhe=cosd(y2) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x2=2*w*costhe*sind(RA/2) x2=x2/15.0 x2=x2+RAmid y2=w*sind(y2) endif call sm_relocate(x2,y2) call sm_expand(0.60) if(x.ge.0) label=x if(x.lt.0) label=abs(abs(x)-24) encode(2,'(i2)',label1) label call sm_putlabel(4,label1) call sm_label('^h') if(x.ge.0) label=(x-label)*60 if(x.lt.0) label=(abs(abs(x)-24)-label)*60 encode(2,'(i2)',label1) label if(label.gt.0) then call sm_label(label1) call sm_label('^m') endif call sm_expand(1.0) endif 17 continue 149 sizedec=DECmax-DECmin sizera=(xmax-xmin)*15.0 if(sizedec.le.1) ddec=0.2 if(sizedec.gt.1.and.sizedec.le.5) ddec=0.5 if(sizedec.gt.5.and.sizedec.le.10) ddec=1.0 if(sizedec.gt.10.and.sizedec.le.20) ddec=2.0 if(sizedec.gt.20.and.sizedec.le.50) ddec=5.0 if(sizedec.gt.50.and.sizedec.le.100) ddec=10.0 if(sizedec.gt.100) ddec=30.0 do 20 x=-90,90,ddec if(x.ge.DECmin.and.x.le.DECmax) then x1=RAmax x2=RAmin y1=x y2=x if(proj.eq.'a') then RA=x1-RAmid RA=RA*15.0 costhe=cosd(y1) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x1=2*w*costhe*sind(RA/2) x1=x1/15.0 x1=x1+RAmid y1=w*sind(y1) RA=x2-RAmid RA=RA*15.0 costhe=cosd(y2) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x2=2*w*costhe*sind(RA/2) x2=x2/15.0 x2=x2+RAmid y2=w*sind(y2) endif call sm_relocate(x1,y1) dx=0.2 if(proj.eq.'a') then do 21 j=1,10 x1=RAmax+((j*dx)*sizera)/170.0/15.0 if(otation.eq.'l') x1=RAmax+((j*dx)*sizera)/227.0/15.0 y1=x RA=x1-RAmid RA=RA*15.0 if(proj.eq.'a') then costhe=cosd(y1) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x1=2*w*costhe*sind(RA/2) x1=x1/15.0 x1=x1+RAmid y1=w*sind(y1) endif call sm_draw(x1,y1) 21 continue endif if(proj.eq.'l') then if(otation.eq.'p') call sm_draw(x1+2.0*sizera/170.0/15.0,y1) if(otation.eq.'l') call sm_draw(x1+2.0*sizera/227.0/15.0,y1) endif dx=0.2 call sm_relocate(x2,y2) if(proj.eq.'a') then do 22 j=1,10 x2=RAmin-((j*dx)*sizera)/170.0/15.0 if(otation.eq.'l') x2=RAmin-((j*dx)*sizera)/227.0/15.0 y2=x RA=x2-RAmid RA=RA*15.0 if(proj.eq.'a') then costhe=cosd(y2) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(RA/2))) x2=2*w*costhe*sind(RA/2) x2=x2/15.0 x2=x2+RAmid y2=w*sind(y2) endif call sm_draw(x2,y2) 22 continue endif if(proj.eq.'l') then if(otation.eq.'p') call sm_draw(x2-2.0*sizera/170.0/15.0,y2) if(otation.eq.'l') call sm_draw(x2-2.0*sizera/227.0/15.0,y2) endif call sm_relocate(x2+(8.0*sizera)/170.0/15.0,y2) if(otation.eq.'l') - call sm_relocate(x2+(8.0*sizera)/227.0/15.0,y2) call sm_expand(0.60) if(ddec.ge.1) then label=x encode(3,'(i3)',label2) label call sm_putlabel(5,label2) call sm_label('^\\circ') endif if(ddec.lt.1) then encode(5,'(f5.1)',label3) x call sm_putlabel(5,label3) call sm_label('^\\circ') endif call sm_expand(1.0) endif 20 continue CCCCCCCCCCCCCCCCC C C C RA/DEC Labels C C C CCCCCCCCCCCCCCCCC if(xmax.gt.xmin) sizera=(xmax-xmin)*15.0 if(xmax.lt.xmin) sizera=-((xmax+24)-xmin)*15.0 sizedec=ymax-ymin RA=(RAmin+RAmax)/2.0 DEC=(DECmin+DECmax)/2.0 if(proj.eq.'a') then x1=RA x1=x1-RAmid x1=x1*15.0 y1=DECmin x2=RAmin x2=x2-RAmid x2=x2*15.0 y2=DEC costhe=cosd(y1) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(x1/2))) x1=2*w*costhe*sind(x1/2) x1=x1/15.0 x1=x1+RAmid y1=w*sind(y1) costhe=cosd(y2) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(x2/2))) x2=2*w*costhe*sind(x2/2) x2=x2/15.0 x2=x2+RAmid y2=w*sind(y2) endif if(proj.eq.'l') then x1=RA x2=RAmin y1=DECmin y2=DEC endif if(abs(RAmax-RAmin).ge.24) goto 170 call sm_relocate(x1,y1-9.0*sizedec/170.0) call sm_putlabel(5,'RA') 170 if(abs(DECmax-DECmin).ge.160) goto 171 if(otation.eq.'p') call sm_relocate(x2+17.0*sizera/170.0/15.0,y2) if(otation.eq.'l') call sm_relocate(x2+17.0*sizera/227.0/15.0,y2) call sm_putlabel(5,'DEC') 171 sizera=sizera/15.0 call sm_relocate(RAmin,DECmin-17.0*sizedec/170.0) if(epoch.eq.2000) call sm_label('Coordinates are J2000.0') if(epoch.eq.1950) call sm_label('Coordinates are B1950.0') if(epoch.ne.1950.and.epoch.ne.2000) then encode(6,'(f6.1)',label3) epoch call sm_relocate(RAmin,DECmin-17.0*sizedec/170.0) call sm_label('Coordinates are') call sm_relocate(RAmin-34.0*sizera/170.0,DECmin-17.0*sizedec/170. -0) call sm_label(label3) endif sizera=sizera*15.0 CCCCCCCCCC C C C Legend C C C CCCCCCCCCC if(legend.eq.'n') goto 150 k=50 j=36 call sm_ctype('default') call sm_ptype(200.0,1) do 25 x=1,7 if(x.eq.1) call sm_expand(minscale+2.5) if(x.eq.2) call sm_expand(minscale+1.5) if(x.eq.3) call sm_expand(minscale+1.0) if(x.eq.4) call sm_expand(minscale+0.8) if(x.eq.5) call sm_expand(minscale+0.5) if(x.eq.6) call sm_expand(minscale+0.25) if(x.eq.7) call sm_expand(minscale) call sm_points(xmin+14.0*sizera/170.0/15.0, - lmin+(2.0+j)*sizedec/170.0,1) call sm_expand(0.4) call sm_relocate(xmin+12.0*sizera/170.0/15.0, - lmin+(2.0+j)*sizedec/170.0,1) if(x.eq.1) call sm_label(' mag < 10') if(x.eq.2) call sm_label(' 10 \\le mag < 12') if(x.eq.3) call sm_label(' 12 \\le mag < 13') if(x.eq.4) call sm_label(' 13 \\le mag < 14') if(x.eq.5) call sm_label(' 14 \\le mag < 15') if(x.eq.6) call sm_label(' 15 \\le mag < 16') if(x.eq.7) call sm_label(' mag \\ge 16') j=j-6 k=k-7 25 continue call sm_expand(1.0) j=0 k=10 do 30 x=1,5 if(x.eq.1) call sm_ptype(203.0,1) if(x.eq.2) call sm_ptype(200.0,1) if(x.eq.3) call sm_ptype(81.0,1) if(x.eq.4) call sm_ptype(41.0,1) if(x.eq.5) call sm_ptype(30.0,1) call sm_expand(0.7) call sm_points(xmin+14.0*sizera/170.0/15.0, - lmax-(1.0+j)*sizedec/170.0,1) call sm_expand(0.4) call sm_relocate(xmin+12.0*sizera/170.0/15.0, - lmax-(1.0+j)*sizedec/170.0,1) if(x.eq.1) call sm_label('Type = -7 to -3') if(x.eq.2) call sm_label('Type = -2 to 0') if(x.eq.3) call sm_label('Type = 1 to 5') if(x.eq.4) call sm_label('Type = 6 to 10') if(x.eq.5) call sm_label('Type = other or blank') j=j+5 k=k+5 30 continue call sm_expand(1.0) lmid=(lmax+lmin)/2.0 if(proj.eq.'a') then x1=xmin+15.0*sizera/170.0/15.0 x1=x1-RAmid x1=x1*15.0 y1=lmid costhe=cosd(y1) w=sqrt((2*(180/PI)*(180/PI))/(1+costhe*cosd(x1/2))) x1=2*w*costhe*sind(x1/2) x1=x1/15.0 x1=x1+RAmid y1=w*sind(y1) endif if(ans2.eq.'y') then k=55 j=15 call sm_ptype(200.0,1) do 31 x=1,colors call sm_ptype(203.0,1) call sm_ctype(c1(x+1)) call sm_points(xmin+14.0*sizera/170.0/15.0, - lmid+j*sizedec/170.0,1) call sm_ctype('default') legendv=v(x) encode(6,'(i6)',label3) legendv call sm_expand(0.4) call sm_label(label3) if(x.eq.1) call sm_label(' \\le V \\le ') if(x.ne.1) call sm_label(' < V \\le ') legendv=v(x+1) encode(6,'(i6)',label3) legendv call sm_label(label3) call sm_expand(1.0) j=j-5 k=k+5 31 continue endif call sm_expand(1.0) 150 if(iterm.eq.2.or.iterm.eq.3.or.iterm.eq.4.or.iterm.eq.5) - call sm_hardcopy call sm_alpha print *,' ' if(iterm.eq.1) print *,'hit return to exit' call sm_redraw(0) if(iterm.eq.5) then print *,' "map.ps" has been created' print *,' This file will be written over if not renamed' endif print *," " print *," " 999 stop end