c STAELPLT. UNIX version of STATNELEV/ELEVMAP. 8-95. D. PLOUFF. c Plots station name and elevation on topographic map base. c CALCOMP plot calls, but PostScript conversion substitutes. c Defaults to 1:24K scale. Map/file names follow Alaska convention but c 1x2-degree sheet basis. Can do Costa Rico 1:50,000 15' maps too. c Allows for paper stretch with requested scale. Emphasize minimum c responses. Option for station names at angles for E-W lines, e.g.. c Options for CB2 or locations only. No large print file now. c PLOUFF GRAVITY DATA format ONLY. METRIC ELEVATIONS PERMITTED. character den*8,dent(6500)*8,dat*9,ans*1 integer INDX(6500),IND(6500),INDEX(6500) dimension ELEV(6500),R(6500),S(6500),STRP(6) common /SORTBK/ XX(6500),KEY(6500) c Save likely repeated parameters except "ifmet" for metric elevs. common /save/ iscale,intv,ifsec,mapn,mapw data LTSM,LTSS,LTND,LTNM,LTNS,LNED,LNEM,LNES,LNWD,LNWM,LNWS, 1 ifutm,intv,ifsec,ifmet,nsta,nmap,ifopen,nosta/19*0/ data dat/' '/,scal,ang45,zero1,zero2/4*0.0/ c "date" requires VAX VMS option -lV77 in UNIX compilation call date (dat) c General prompt. Data input file/parameters. Loop individual maps. call promg (ifmet,ifutm,dat) if (ifutm .gt. 1) stop C PLOUFF format--test input file for format errors. CB2 to 0.1 mGal. 1 if (ifmet .ne. 8) then read (9,900,end=3,err=2) DEN,TDEG,TMIN,GDEG,GMIN,EL 900 format (bz,A8,1X,F2.0,F4.2,1X,F3.0,F4.2,F6.1) else read (9,908,end=3,err=2) DEN,TDEG,TMIN,GDEG,GMIN,EL 908 format (bz,A8,1X,F2.0,F4.2,1X,F3.0,F4.2,46x,F6.1) end if nsta=nsta+1 go to 1 2 nsta=nsta+1 print 601, nsta write (16,601) nsta 601 format (' **STOP. Bad format on input line',I5) go to 99 3 write (16,602) nsta print 602, nsta 602 format (I5,' stations in input file.',/) rewind 9 c Re-read input file later for each map. Parameter error test. 5 ltsd=99 C Read prompted parameters for next map to be plotted call promp (LTSD,LTSM,LTSS,LTND,LTNM,LTNS,LNED,LNEM,LNES, 1 LNWD,LNWM,LNWS,IFMET,nmap,scal,ang45,nosta) IF (LTSD .GT. 90) go to 5 C INITIALIZE PLOTTER; MOVE ORIGIN AWAY FROM EDGE. DRAW L-SHAPED MARK C FOR MECHANICAL PEN ORIGIN 0.5" SOUTH AND 0.25" WEST OF SW MAP CORNER. call plots (0,0,0) ifopen=1 call plot (0.25,0.5,-3) c Edges and intervals in seconds. Old version had minute option. KTSM=(60*LTSD+LTSM)*60+LTSS KTNM=(60*LTND+LTNM)*60+LTNS KNEM=(60*LNED+LNEM)*60+LNES KNWM=(60*LNWD+LNWM)*60+LNWS c Number of intervals in each direction KNORT=(KTNM-KTSM)/intv KWEST=(KNWM-KNEM)/intv c Number of tic marks in each direction and total JNORT=KNORT+1 JWEST=KWEST+1 N=JNORT*JWEST c Convert meters to inches at given scale reciprocal SCALE=39.37008/SCAL FTIK =float(intv)/3600.0 DSOUT=float(KTSM)/3600.0 DEAST=float(KNEM)/3600.0 DNORT=float(KTNM)/3600.0 DWEST=float(KNWM)/3600.0 c Meridian & latitude at center of map--orients and maximizes accuracy CEN=(KNEM+KNWM)/7200.0 CLAT=(KTSM+KTNM)/7200.0 c Initialize projection constants for mid-latitude and zero meridian call proj (ifutm,0,clat,0.0,zero1,zero2,0.0) L=0 do 6 J=1,JNORT dpi=DSOUT+FTIK*(J-1)-CLAT do 6 K=1,JWEST L=L+1 c East to west in rows that progress from south to north dli=CEN-(DEAST+FTIK*(K-1)) c Slightly greater time for latitude operations in longitude loop call proj (ifutm,1,dpi,dli,s(L),r(L),scale) 6 continue c All coordinates in inches relative to southwest SMAX=-S(jwest) WMAX=-R(jwest) c Protect against exceeding array limit of 6500 in prompting do 8 J=1,N R(J)=R(J)+WMAX 8 S(J)=S(J)+SMAX JWEST1=JWEST+1 KEAST=N-KWEST YKE=S(KEAST) X1=R(1) Y1=S(1) XKE=R(KEAST) XN=R(N) YN=S(N) C PRINT/PLOT INPUT FILE NAME (not needed) AND DATE. call symbol's c(-0.4,0.04,0.07,'STATION FILE:',90.0,13)(-0.4,1.0,0.07,infile,90.0,65) call symbol (-0.15,0.1,0.07,DAT,90.0,9) if (nosta .eq. 2) go to 4 if (ifmet .eq. 0) call symbol (-0.15,0.75,0.07,'*FEET',90.0,5) if (ifmet .eq. 1) call symbol (-0.15,0.75,0.07,'*METERS',90.0,7) if (ifmet .eq. 8) call symbol (-0.15,0.75,0.07,'0.1-MGAL CB2', 1 90.0,12) C LINES ALONG EDGES AT MECHANICAL PEN ORIGIN. Print at 30 degrees. 4 call symbol (-0.19,-0.47,0.07,'ORIGIN',30.0,6) call plot (-0.15,0.0,3) C X=-0.5 (MECHANICAL ZERO) DID NOT DRAW; TRY -0.496 (-0.248 REL -.25) call plot (-0.248,0.0,2) call plot (-0.248,-0.5,2) call plot (0.,-0.5,2) call plot (0.,-0.2,2) call plot (0.,-0.5,2) C REPEAT POINT call plus (0.1,-0.35,0.16) call plus (0.1,-0.35,0.16) c First, draw tic marks. Start in SW corner call ell(r(jwest),s(jwest)+0.12,r(jwest),s(jwest), 1 r(jwest)+0.12,s(jwest)) c if only one east-west interval on the map if (kwest .lt. 2) go to 25 do 27 J=2,KWEST K=JWEST1-J 27 call tee(r(k)-0.12,s(k),r(k),s(k),r(k),s(k)+0.12,r(k)+0.12,s(k)) 25 call ell (r(1)-0.12,s(1),r(1),s(1),r(1),s(1)+0.12) call geog (LTSD,LTSM,LTSS,X1+0.16,Y1-0.07,IFSEC) c If only one north-south interval on the map if (knort .eq. 1) go to 69 LNORTM=KNORT-1 c Sets direction of movement. Negative is left to right. LIN=-1 do 30 L=1,LNORTM LIN=-LIN IF (LIN .lt. 0) go to 29 c Right to left JB=L*JWEST+1 KE=JB+KWEST call tee (r(jb),s(jb)-0.12,r(jb),s(jb),r(jb)-0.12,s(jb),r(jb), 1 s(jb)+0.12) if (kwest .lt. 2) go to 32 jbp=jb+1 kem=ke-1 do 67 k=jbp,kem 67 call plus (r(k),s(k),0.24) 32 call tee (r(ke),s(ke)-0.12,r(ke),s(ke),r(ke)+0.12,s(ke),r(ke), 1 s(ke)+0.12) go to 30 c Left to right 29 jb=(l+1)*jwest call tee (r(jb),s(jb)-0.12,r(jb),s(jb),r(jb)+0.12,s(jb),r(jb), 1 s(jb)+0.12) jb=jb+1 if (kwest .lt. 2) go to 41 do 68 j=2,kwest k=jb-j 68 call plus (R(K),S(K),0.24) 41 jb=jb-jwest call tee (r(jb),s(jb)-0.12,r(jb),s(jb),r(jb)-0.12,s(jb),r(jb), 1 s(jb)+0.12) 30 continue 69 call ell (xke,yke-0.12,xke,yke,xke-0.12,yke) call geog (LTND,LTNM,LTNS,XKE+0.16,YKE-0.07,IFSEC) call geog (LNED,LNEM,LNES,XKE-0.36,YKE+0.13,IFSEC) KE1=KEAST+1 if (kwest .lt. 2) go to 72 nm=n-1 do 71 k=ke1,nm 71 call tee(r(k)+0.12,s(k),r(k),s(k),r(k),s(k)-0.12,r(k)-0.12,s(k)) 72 call ell (r(n),s(n)-0.12,r(n),s(n),r(n)+0.12,s(n)) call geog (LNWD,LNWM,LNWS,XN-0.23,YN+0.13,IFSEC) XLAST=XN YLAST=YN M=0 C Initialize unlikely distances. Saves one test in DO-loop. SMX=5000.0 BGX=-50.0 do 7 J=1,NSTA if (ifmet .ne. 8) read (9,900) DEN,TDEG,TMIN,GDEG,GMIN,EL c To feet or 0.1 mGal if (ifmet .eq. 8) read (9,908) DEN,TDEG,TMIN,GDEG,GMIN,EL PLAT=TDEG+TMIN/60.0 IF (PLAT .lt. DSOUT .or. plat .gt. dnort) go to 7 PLON=ABS(GDEG)+GMIN/60.0 IF (plon .lt. deast .or. plon .gt. dwest) go to 7 DP=PLAT-CLAT M=M+1 IF (M .gt. 6500) then print 603 write (16,603) 603 format (' **SKIP MAP WITH MORE THAN 6500 DATA POINTS.',/, 1 ' **ONLY TICKMARKS WERE PLOTTED.') go to 999 end if dl=cen-plon call proj (ifutm,1,dp,dl,yj,xj,scale) xj=xj+wmax yj=yj+smax R(M)=XJ S(M)=YJ dent(m)=den elev(m)=el IF (XJ .GT. BGX) BGX=XJ IF (XJ .LT. SMX) SMX=XJ 7 continue if (m .eq. 0) then print 604 write (16,604) 604 format (' This map has no data points on the plot.',/) nst=0 go to 15 end if INDEX(1)=1 C STRIP WIDTHS FOR OPTIMUM PEN MOVEMENT 2.0 TO 4.5 INCH AT 0.5 INCH STRIP=2.0 do 26 JSTR=1,6 STRP(JSTR)=STRIP STRIP=STRIP+0.5 c 26 DIST(JSTR)=0.0 26 IF (M .eq. 1) go to 33 do 23 JSTR=1,6 STRIP=STRP(JSTR) XBEG=XLAST YBEG=YLAST DS=0.0D0 JN=0 NST=1 C NUMBER OF STRIPS JDIFF=1.0+(BGX+0.01-SMX)/STRIP XMIN= SMX-0.001 KASC= 1 do 35 L=1,JDIFF QMIN=XMIN+STRIP NST=NST+JN C JN-NUMBER OF POINTS IN STRIP. NST-ONE MORE THAN TOTAL OF POINTS THAT C INCLUDE LAST STRIP. KASC-POSITIVE FOR ASCENDING ORDER. JN=0 do 36 J=1,M XJ=R(J) c IF (XJ .GT. XMIN .AND. XJ .LE. QMIN) go to 38 IF (XJ .le. XMIN .or. XJ .gt. QMIN) go to 36 JN=JN+1 INDX(JN)=J XX(JN)=S(J) 36 continue XMIN=QMIN IF (JN-1) 35,79,40 79 KEY(1)=1 go to 37 40 KASC=-KASC call sort (JN,KASC) 37 NSTM=NST-1 do 42 J=1,JN JL=NSTM+J K=KEY(J) LND=INDX(K) IND(JL)=LND XJ=R(LND) YJ=S(LND) DS=DS+SQRT((XJ-XBEG)**2+(YJ-YBEG)**2) XBEG=XJ 42 YBEG=YJ 35 continue NST=NST+JN-1 IF (NST .EQ. M) go to 24 c Path only if a later erroneous change of the algorithm. write (16,605) NST,M print 605, NST,M 605 format(' **STOP. TOTAL POINTS SORTED,',I5,' NOT TOTAL ON MAP,', 1 I5) nst=0 go to 15 c Previously skipped to 34 if number of strips was one rather than 6. C PEN STARTED AT NW CORNER; STOPS AT REPEAT POINT (SW) 24 XJ=XJ-0.2 DS=DS+SQRT(XJ*XJ+YJ*YJ) c DIST(JSTR)=DS; c Initialize if first of the loop of 6 strips. IF (JSTR .EQ. 1) go to 34 c Don't reset indices if this path is longer than shortest to now. IF (DS .GE. DMIN) go to 23 34 DMIN=DS do 39 I=1,M 39 INDEX(I)=IND(I) 23 continue c write (STRP(I),DIST(I),I=1,6); format (' LENGTHS',6(F4.1,F6.0)) 33 do 14 I=1,M C INDEX IN R,S data point array LND=INDEX(I) XJ=R(LND) YJ=S(LND) call plus (XJ,YJ,0.06) EL=ELEV(lnd) c Metric units for elevations if (ifmet .eq. 1) el=0.3048*el c Protect against elevation truncation instead of rounding jel=abs(el)+0.5 if (el .lt. 0.0) jel=-jel el=jel if (nosta .ne. 0) go to 14 c Fonts (except Courier) change consistency of alinement rel "+". if (ang45 .eq. 0.0) call symbol (xj-0.24,yj-0.12,0.07,dent(lnd), 1 0.0,8) if (ang45 .ne. 0.0) call symbol (xj-0.07,yj-0.12,0.07,dent(lnd), 1 ang45,8) 14 if (nosta .ne. 2) call number (XJ+0.07,YJ+0.04,0.07,EL,90.0,-1) c Plot number of stations plotted fm=m call symbol (-0.15,1.45,0.07,'STATIONS=',90.0,9) call number (-0.15,2.03,0.07,fm,90.0,-1) C REPEAT POINT 15 call symbol (0.2,-0.38,0.07,'REPEAT ',0.0,12) c call square (0.1,-0.35,0.04) call circle (0.1,-0.35,0.05) if (ifutm .eq. 1) then call symbol (1.5,-0.38,0.07,'UTM 1:',0.0,6) call number (1.87,-0.38,0.07,scal,0.0,-1) else call symbol (1.5,-0.38,0.07,'POLYCONIC 1:',0.0,12) call number (2.26,-0.38,0.07,scal,0.0,-1) end if write (16,609) M print 609, M 609 format (i5,' stations on map.',/) C Terminate this plot, rewind data input file to get new subset. 999 call plot (0.0,0.0,999) close (7) ifopen=0 rewind (9) print 101 101 format (' Do you want to plot more maps?') read 501, ans 501 format (a1) if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 5 99 print 699, nmap write (16,699) nmap 699 format (i5,' map files were written.') close (9) close (16) if (ifopen .ne. 0) close (7) stop end c ----------------------------------------------------------- subroutine geog (JD,JM,JS,X,Y,IFSEC) C TO PLOT DEGREES-MINUTES-SECONDS COORDINATE. ORIGIN AT C LOWER LEFT CORNER OF DEGREE VALUE. NUMBER WIDTH 0.14 INCH. Z=JD call number (X,Y,0.14,Z,0.0,-1) XX=X+0.35 IF (JD .GT. 99) XX=XX+0.12 c call square (XX-0.08,Y+0.14,0.04) call circle (XX-0.10,Y+0.13,0.05) Z=JM XNEW=XX+0.34 IF (JM .LT. 10) XX=XX+0.12 call number (XX,Y,0.14,Z,0.0,-1) YS=Y+0.12 call plot (XNEW-0.07,YS,3) call plot (XNEW-0.06,YS+0.04,2) call plot (XNEW-0.07,YS,2) c IFSEC is a flag for no seconds in tickmark interval IF (IFSEC .eq. 0) return Z=JS XX=XNEW+0.34 IF (JS .LT. 10) XNEW=XNEW+0.12 call number (XNEW,Y,0.14,Z,0.0,-1) call plot (XX-0.08,YS,3) call plot (XX-0.06,YS+0.04,2) call plot (XX-0.03,YS+0.04,3) call plot (XX-0.05,YS,2) return end c ------------------------------------------------------------------- subroutine sort (NO,IND) C SORT ARRANGES THE ELEMENTS OF X IN ASCENDING OR DESCENDING ORDER AND C CONSTRUCTS AN ARRAY KEY OF SUBSCRIPTS OF X ARGUMENTS. 'SHELL SORT' c -IBM VERSION. NO =THE NUMBER OF ELEMENTS IN THE X ARRAY C IND=1 IF ASCENDING ORDER; =-1 IF DESCENDING ORDER COMMON /SORTBK/ X(6500),KEY(6500) do 1 I=1,NO 1 KEY(I)=I IF(IND.EQ.1) go to 5 do 3 I=1,NO 3 X(I)=-X(I) 5 MO=NO 2 IF(MO .GT. 15) go to 23 IF(MO .LT. 2) go to 9 MO=2*(MO/4)+1 go to 24 23 MO=2*(MO/8)+1 24 KO=NO-MO JO=1 25 I=JO 26 TEMP=X(I) HEMP=X(I+MO) IF (HEMP .GE. TEMP) go to 28 X(I)=HEMP X(I+MO)=TEMP KEMP=KEY(I) KEY(I)=KEY(I+MO) KEY(I+MO)=KEMP I=I-MO IF(I .GT. 0) go to 26 28 JO=JO+1 IF(JO-KO)25,25,2 9 IF(IND.EQ.1) return do 10 I=1,NO 10 X(I)=-X(I) return end c -------------------------------------------------------- subroutine plus (X,Y,H3) C DRAW + SIGN REFERRED TO 4/6 OF INCREMENT HEIGHT H=H3/3.0 XM=X-H XP=X+H YM=Y-H YP=Y+H call plot (XM,Y,3) call plot (XP,Y,2) call plot (XM,Y,2) call plot (X,YM,3) call plot (X,YP,2) call plot (X,YM,2) return end c ------------------------------------------------------------- c subroutine square (X,Y,H2) C DRAW A SQUARE for repeat point and degree symbols c H=0.5*H2; plot (X-H,Y-H,3); plot (X-H,Y+H,2); plot (X+H,Y+H,2) c plot (X+H,Y-H,2); call plot (X-H,Y-H,2); return; end c ------------------------------------------------------------- subroutine promg (ifmet,ifutm,dat) c General prompting and file opening. Specific maps in PROMP character infile*65,ans*1,dat*9 logical exists print 100 100 format (' STAELPLT. Plouff, 4-98. Plot station names and ', 1 'elevations or anomaly 2 on',/,' topo maps. Input data are ', 2 'in plouff format. Default maps 7-1/2-minutes',/,' at ', 3 'scale of 1:24,000. Map file names based on 15-minute ', 4 'naming convention',/,' (A-D south to north within a 1- by ', 5 '2-degree quadrangle and 1-8 east to the',/,' west and NE-', 6 'NW-SE-SW for 7.5 minutes) are suggested.',/,' Responses of ', 7 'y, Y, or a carriage return signify yes.') inquire (file='STAELPLT.PNT',exist=exists) if (exists) then print 101 101 format (' A file for recording this session, STAELPLT.PNT, ', 1 'already exists.',/,' Do you want to STOP to save it?') read 501, ans 501 format (a1) if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') then ifutm=9 return end if end if print 102 102 format (' TYPE the name of your input file (plouff format):') read 565, infile 565 format (a65) inquire (file=infile,exist=exists) if (.not. exists) then print 103 103 format (' That file does not exist. Try once more.') print 102 read 565, infile inquire (file=infile,exist=exists) if (.not. exists) then print 104 104 format (' **STOP. File again was not found.') ifutm=9 return end if end if print 105 105 format (' Do you want to plot station elevations (not CB2)?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 1 ifmet=8 print 106 106 format (' Anomaly 2 to 0.1 mGal (cols 70-75) will be plotted ', 1 'on all maps.') 1 print 107 107 format (' Map plots must be with the transverse mercator ', 1 'or polyconic projection.',/,' Is your projection the common', 2 ' ''UTM'' (transverse mercator)?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') ifutm=1 open (9,file=infile,form='formatted',status='old',blank='zero') open (16,file='STAELPLT.PNT',form='formatted',status='unknown') write (16,160) dat,infile 160 format (1x,a9,' Program STAELPLT.'/,' Data file: ',a65) if (ifutm .eq. 0) write (16,161) 161 format (' All maps have the polyconic projection.') if (ifutm .eq. 1) write (16,162) 162 format (' All maps have the transverse mercator (''UTM'') ', 1 'projection.') if (ifmet .eq. 8) write (16,106) return end c ------------------------------------------------------------------ subroutine promp (ltsd,ltsm,ltss,ltnd,ltnm,ltns,lned,lnem, 1 lnes,lnwd,lnwm,lnws,ifmet,nmap,scal,ang45,nosta) character*1 ans,angles(8) c Station name at 10-degree intervals for "ang45" angle & -45 degrees data angles/'1','2','3','4','5','6','7','8'/ c Save likely repeated parameters except "ifmet" for metric elevs. common /save/ iscale,intv,ifsec,mapn,mapw if (nmap .ne. 0) then print 100 100 format (' Unless there was a mistake on the last attempt,',/, 1 ' Do you want the same scale, tickmark interval, and map ', 2 'size as the last map?') read 501, ans 501 format (a1) if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 10 end if print 101 101 format (' Is this a conventional 1:24,000 map?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') then iscale=24000 intm=2 ifsec=30 intv=150 mapn=30+7*60 mapw=30+7*60 go to 10 else print 102 102 format (' Is this a conventional 1:62,500 map?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') then iscale=62500 intm=5 ifsec=0 intv=300 mapn=15*60 mapw=15*60 go to 10 end if end if 3 print 103 103 format (' TYPE the reciprocal of the map scale (INTEGER with', 1 'out commas) (e.g. 50000):') read (5,*,err=3) iscale 4 print 104 104 format (' TYPE the tickmark interval (INTEGER minutes, sec', 1 'onds) (e.g. 2 30):') read (5,*,err=4) intm,ifsec c Later, IFSEC will be flag for suppression of seconds label intv=ifsec+60*intm if ((3600-intv*(3600/intv)) .ne. 0) then print 105 105 format (' Try again. One degree must be an INTEGER multi', 1 'ple of your tickmark interval.') go to 4 end if 6 print 106 106 format (' TYPE the east-west map width (INTEGER minutes, sec', 1 'onds) (e.g. 7 30):') read (5,*,err=6) imw,isw mapw=60*imw+isw if ((mapw-intv*(mapw/intv)) .ne. 0) then print 107 107 format (' Try again. Your map width is not an INTEGER multi', 1 'ple of your tickmark interval.') go to 4 end if mapn=mapw print 108 108 format (' Is the north-south map width the same as the east-', 1 'west?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 9 8 print 109 109 format (' TYPE the north-south map width (INTEGER minutes, sec', 1 'onds) (e.g. 10 0):') read (5,*,err=8) imw,isw mapn=60*imw+isw if ((mapn-intv*(mapn/intv)) .ne. 0) then print 107 go to 8 end if 9 n=(1+mapw/intv)*(1+mapn/intv) if (n .le. 6500) go to 10 print 117, n 117 format (i10,' exceeds tickmark limit of 6500. Try again.') go to 4 10 print 110 110 format (' The southwest corner must be typed with INTEGERS (de', 1 'grees/minutes/seconds)',/,' separated by blank spaces.',/, 2 ' TYPE 999 for degrees to correct an earlier mistake.') 11 print 111 111 format (' TYPE the south edge of map (3 INTEGERS):') read (5,*,err=11) ltsd,ltsm,ltss ltsd=iabs(ltsd) if (ltsd .gt. 89) go to 3 call findedge (1,mapn,ltsd,ltsm,ltss,ltnd,ltnm,ltns) if (ltsd .eq. 999) go to 11 12 print 112 112 format (' TYPE west edge of map (3 INTEGERS):') read (5,*,err=12) lnwd,lnwm,lnws lnwd=iabs(lnwd) if (lnwd .gt. 360) go to 11 call findedge (-1,mapw,lnwd,lnwm,lnws,lned,lnem,lnes) if (lnwd .eq. 999) go to 12 scal=float(iscale) iangle=0 print 114 114 format (' A response to the following question provides severa', 1 'l functions:',/,' [y] default--horizontal; numbers [1] to ', 2 '[8]--plot at 10 to 80 degrees SE;',/,' [-]--do not plot sta', 3 'tion names; [=]--only plot locations;',/,' [others]--plot ', 4 'at 45 degrees SE.',//,' Horizontal station names OKAY?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') then ang45=0.0 else nosta=0 if (ans .eq. '-' .or. ans .eq.'=') then if (ans .eq. '-') then nosta=1 print 115 write (16,115) 115 format (' Station names will not be plotted.') else nosta=2 print 116 write (16,116) 116 format (' Only locations will be plotted.') end if go to 17 end if do 15 i=1,8 if (angles(i) .ne. ans) go to 15 iangle=10*i go to 16 15 continue c Default for E-W lines; answer "n", e.g. No print for horizontal. iangle=45 16 write (16,601) iangle print 601, iangle 601 format (' Station names are plotted at',i3,' degrees.') 17 ang45=-float(iangle) end if c Bypass if CB2 is plotted or only locations are plotted. if (ifmet .eq. 8 .or. nosta .eq. 2) go to 13 ifmet=1 print 113 113 format (' Are elevations on this map expressed in feet (no=', 1 'meters)?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') ifmet=0 c OPEN FILE 7 BEFORE write--HERE RATHER THAN IN call plot(0,0,999) 13 call openmap (ltsm,ltss,lnwd,lnwm,lnws) c No way to stop here. write (16,600) ltsd,ltsm,ltss,ltnd,ltnm,ltns,lned,lnem,lnes, 1 lnwd,lnwm,lnws,iscale,intm,ifsec 600 format (' Map extends from',i3,'D',i3,'M',i3,'S to',i3,'D',i3, 1 'M',i3,'S and',i4,'D',i3,'M',i3,'S to',i4,'D',i3,'M',i3,'S',/, 2 ' at a scale of 1:',i6,' and a tickmark interval of',i3,'M', 3 i3,'S') nmap=nmap+1 return end c ------------------------------------------------------------------ subroutine openmap (ltsm,ltss,lnwd,lnwm,lnws) c Open a map file, suggesting a default map name, here rather than c in subroutine for call plots (0.,0.,0.), as in other applications. character mapfile*20,temfile*20,blank*20,ans*1,tem(20)*1, 1 blnk(20)*1,latcode(4)*1,loncode(8)*1 equivalence (temfile,tem(1)),(blank,blnk(1)) logical exists common /save/ iscale,intv,ifsec,mapn,mapw data blnk,tem/40*' '/,latcode/'A','B','C','D'/ data loncode/'1','2','3','4','5','6','7','8'/ c Allow 5% scalings for map type definitions 24K/25K, 50K/62.5K/63360 if ((iscale .gt. 22800 .and. iscale .lt. 26250) .or. (iscale .gt. 1 47500 .and. iscale .lt. 66500)) then temfile=blank c Latitude minutes (0,7+), (15,22+), (30,37+), or (45,52+) lt=1+ltsm/15 tem(1)=latcode(lt) c Measure from east edge of 2-degree sheet. Adjust east 7 1/2s & sheet. ln=((lnwd-(2*(lnwd/2)))*60+lnwm)/15 if ((lnwm-(15*(lnwm/15))) .ne. 0) ln=ln+1 if (ln .eq. 0) ln=8 tem(2)=loncode(ln) if (iscale .gt. 30000) then tem(3)='.' else if (ltss .eq. 30) tem(3)='N' if (ltss .eq. 0) tem(3)='S' if (lnws .eq. 30) tem(4)='E' if (lnws .eq. 0) tem(4)='W' c SW corner not for a standard 7.5-minute quad if (tem(4) .eq. ' ' .or. tem(3) .eq. ' ' ) go to 9 tem(5)='.' end if mapfile=temfile print 100, mapfile 100 format (' An appropriate plot file name for your map is ',a20, 1 /,' Is this name acceptable?') read 501, ans 501 format (a1) if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 10 end if 9 print 109 109 format (' TYPE a plot file name for this map (<21 cols):') read 520, mapfile 520 format (a20) 10 inquire (file=mapfile,exist=exists) if (exists) then print 110 110 format (' That plot file already exists.',/,' Do you want ', 1 'to overwrite that plot file?') read 501, ans if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') go to 99 go to 9 end if 99 open (7,file=mapfile,form='formatted',status='unknown') write (16,600) mapfile 600 format (' Map file name is: ',a20) return end c ------------------------------------------------------------------ subroutine findedge (isign,map,ld,lm,ls,kd,km,ks) c Add or subtract one map width in geographic coordinates ks=ls+60*lm+3600*ld itest=ks-map*(ks/map) if (itest .eq. 0) go to 2 print 500 500 format (' Map is not an integral number of map-widths from ', 1 'the Equator or Greenwich.') c ld=999; return 8-14-95 WARNING, NOT STOP RESTRICTION c isign--increase toward north; decrease toward east from SW corner. 2 ks=ks+isign*map km=ks/60 ks=ks-60*km kd=km/60 km=km-kd*60 return end c --------------------------------------------------------------- subroutine ell(x1,y1,x2,y2,x3,y3) c Draw L-shaped line call plot (x1,y1,3) call plot (x2,y2,2) call plot (x3,y3,2) return end c ------------------------------------------------------------- subroutine tee(x1,y1,x2,y2,x3,y3,x4,y4) c Draw T-shaped line call plot (x1,y1,3) call plot (x2,y2,2) call plot (x3,y3,2) call plot (x2,y2,2) call plot (x4,y4,2) return end c ------------------------------------------------------------ subroutine proj (ifutm,ifcons,DPI,DLI,DP,DL,scale) common /params/ age,agef,cp,sp,fm1,sc228,sc5 c UTM and Polyconic conversions (Plouff, 1968, p. C175)--PP600C DPR=1.745329E-2*DPI DLB=1.745329E-2*DLI if (ifcons .eq. 1) then C INPUT DL IS CENTRAL MERIDIAN OF MAP MINUS LONGITUDE IN DEGREES. C INPUT DP IS LATITUDE MINUS CENTRAL LATITUDE OF MAP, IN DEGREES. C OUTPUT Y-COORDINATE IS DP. X IS DL. BOTH ARE IN meters. if (ifutm .ne. 0) then c Routine UTM conversion GL=(AGE+DPR*AGEF)*DLB DL=scale*GL*(CP-DPR*SP) DP=scale*(DPI*(FM1+DPR*SC228)+DLB*GL*SC5) else c Routine Polyconic conversion DL=scale*DLB*(AGEF-DPR*SP) DP=scale*(DPI*FM1+DLB*DLB*AGE) end if else c Input for constants is latitude at center of map--CLAT=DPI CP=COS(DPR) SP=SIN(DPR) CP2=CP*CP if (ifutm .ne. 0) then c Assign UTM constants SP2=1.0-CP2 SCP=SP*CP SC5=0.5*SCP F1141=1141.70-9.60*CP2 FM1=(111699.34-CP2*F1141)*0.9996 SC228=F1141*SCP*0.9996 GE=1.0+0.5*6.768658E-3*SP2*(1.0+0.75*6.768658E-3*SP2) AGE=6378206.4*GE*0.9996 AGEF=6378206.4*6.768658E-3*SCP*0.9996 else c Assign Polyconic constants SCP=0.5*SP*CP FM1=111699.34+CP2*(9.60*CP2-1141.70) AGE=6378206.4*SCP AGEF=(6.768658E-3*SCP*SP+CP)*6378206.4 SP=SP*6378206.4 end if end if return end c ---------------------------------------------------------------- subroutine plots (idum,jdum,ldev) c D. Plouff 12-94. Subroutines to convert CALCOMP plot calls to c PostScript for STAELPLT program. Open plot file. c Plotter initialization routine - by MAIN before other plot calls. c idum, jdum, ldev = dummy variables for compatibility with caller. c Test and open plot file with name and unit 7 outside. c Set up the header required to identify following as PS commands. character*1 ans write (7,700) 700 format ('%!PS-Adobe') print 500 500 format (' PostScript-2 printers require a BoundingBox specific', 1 'ation. The default',/,' size is 8.5 by 11 inches. Others ', 2 'will be 11 by 17 inches (-tabloid). The',/,' BoundingBox ', 3 'line number 2 in the output file can be edited to yield a ', 4 'size',/,' larger than 11 by 17 (1 inch = 72 points).') print 501 501 format (' Is the size east 8.5 by north 11 inches okay?') read 502, ans 502 format (a1) if (ans .eq. 'y' .or. ans .eq. 'Y' .or. ans .eq. ' ') then write (7,701) 701 format ('%%BoundingBox 50 50 600 780') else write (7,702) 702 format ('%%BoundingBox 50 50 780 1200') end if c Bounding box to save paper on some plotters "cc" did not work "<" cc write (7,701) cc 701 format ('%%EndComments',/,'%%BeginDefaults') c Set up definitions for other routines: M, change to a new location. c S, print previous character string in parentheses. write (7,703) 703 format ('/M {moveto} def /S {show} def',/, 1 '/LTSM {lineto currentpoint stroke M} def') c After definitions, before scale change. [2448 3168] for Davinci? c One initial Davinci translation is "500 3000" not "0 0"--Bob Morin cc write (7,703) cc 703 format ('%%BeginSetup',/,'<<',/,' /InputAttributes <<',/, cc 1 ' 0 << /Pagesize [612 792] >>',/,' >>',/,'<< setpage', cc 2 'device',//,'<< /PageSize [612 792] >> setpagedevice',/, cc 3 '%%EndSetup',/,'0 0 translate') c Units inches, origin near lower left corner. write (7,704) 704 format ('72 72 scale') c Set line width 300 dpi laser write (7,705) 705 format ('0.004 setlinewidth') c Extra origin shift for page plots write (7,706) 706 format ('0.3 0.3 translate') return end c --------------------------------------------------------------- subroutine plot (x,y,ipen) c CalComp calls. x,y coordinates in inches from the origin. c ipen=2, move/draw with pen down; =3, move with pen up; =999, c move with pen up, terminate plot, close plot file if (ipen .eq. -3) then c ipen=-3 for initial origin change write (7,1) x,y 1 format (2f8.4,' translate') return end if c Pen down - draw if (ipen .eq. 2) then write (7,2) x,y 2 format (2f8.4,' LTSM') c Pen up - move if ipen=3 (start a line) or 999 by default else write (7,3) x,y 3 format (2f8.4,' M') end if if (ipen .eq. 999) then c Print this only page and close plot file write (7,999) 999 format ('showpage % end of plot') close (7) end if return end c --------------------------------------------------------------- subroutine number (x,y,hgt,fpn,angle,ndec) c Numbers at angles=0 and 90 degrees (ndec=-1) only. c Heights 0.06 and 0.12 (tickmark labels only) inch. 0.07--0.113 if (hgt .gt. 0.10) write (7,200) x,y 200 format ('/Helvetica findfont 0.194 scalefont setfont',/, 1 2f8.3,' M') if (hgt .lt. 0.09) write (7,201) x,y 201 format ('/Helvetica findfont 0.097 scalefont setfont',/, 1 2f8.3,' M') if (angle .gt. 89.0) write (7,10) 10 format (' 90.0 currentpoint gsave translate rotate') c Find number of characters. Negative numbers need space for minus. if (fpn .gt. -0.1) iv=fpn+0.1 if (fpn .lt. -0.9) iv=fpn-0.1 if (iv .gt. -1 .and. iv .lt. 10) then write (7,301) iv 301 format ('(',i1,') S') go to 2 end if if (iv .gt. -10 .and. iv .lt. 100) then write (7,302) iv 302 format ('(',i2,') S') go to 2 end if if (iv .gt. -10 .and. iv .lt. 1000) then write (7,303) iv 303 format ('(',i3,') S') go to 2 end if if (iv .gt. -1000 .and. iv .lt. 10000) then write (7,304) iv 304 format ('(',i4,') S') go to 2 end if if (iv .gt. -10000 .and. iv .lt. 100000) then write (7,305) iv 305 format ('(',i5,') S') go to 2 end if if (iv .gt. -100000 .and. iv .lt. 1000000) then write (7,306) iv 306 format ('(',i6,') S') go to 2 end if print 100 100 format (' Number=',i9,' was not plotted.') 2 if (angle .gt. 89.0) write (7,19) 19 format ('grestore') return end c -------------------------------------------------------------------- subroutine symbol (x,y,ht,name,ang,nch) c Set string to max length of 12 (REPEAT POINT;POLYCONIC 1:) character name*12,temp*12,char(12)*1 equivalence (temp,char(1)) c Heights 0.06 inch for this application despite 0.07 call. MAPSCALE c shows slight difference in Helvetica NUMBER and uppercase heights. write (7,100) 100 format ('/Helvetica findfont 0.09 scalefont setfont') c If left-adjusted and padded with blanks, all may print as a12, but c need to replace unprotected trailing characters with blank spaces. temp=name if (nch .lt. 12) then nchp=nch+1 do 1 i=nchp,12 1 char(i)=' ' end if if (ang .lt. 0.1 .and. ang .gt. -0.1) then c Horizontal orientation for station ID and "REPEAT[ POINT]". write (7,102) x,y,temp 102 format (2f8.3,' M (',a12,') S') return end if c Rotated. Check for offset Courier/Helvetica Station -45 degrees. cp yp=y+0.07 yp=y write (7,103) x,yp,ang,temp 103 format (2f8.3,' M',/,f5.1,' currentpoint gsave translate ', 1 'rotate',/,'(',a12,') S',/,'grestore') return end c ---------------------------------------------------------------- subroutine circle (x,y,h2) C Draw a circle for repeat point and degree symbols. 4-98, Plouff. h=0.5*h2 write (7,700) x,y,h 700 format ('newpath ',2f8.4,f6.3,' 0 360 arc stroke') return end