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