c--------------------------------------------------------------------- c c copyright (c) 1986 by c the Bone and Calcium Research Laboratory c c/o Metabolic Unit c Repatriation General Hospital c Concord NSW 2139 c Australia c c This software is furnished under license and may be used c and copied only in accordance with the terms of such c license and only with the inclusion of the above copyright c notice. This software or any other copies thereof may be c provided or otherwise made available to any other person. c No title to and ownership of the software is hereby c transferred. c c Support in the implementation or advice in the usage of c the software cannot be offered by BCRL. Criticism and c comments are nevertheless welcomed. c c The BCRL assumes no responsibility for the use, reliability c or any damage of any nature which may result from the use c of this software. c c--------------------------------------------------------------------- c--------------------------------------------------------------------- c c NAME c aaxis.for c c USAGE c call axxis(xtic,ytic) c c ABSTRACT c Subroutine to draw axis for GRAPH. It's differs c from AXIS in that a) starts axis at the minimum points c rather than the minimum plot window points, and b) on the c plotter offsets 1.4 from the plot area. c c INPUT c xtic - x axis tic interval c ytic - y axis tic interval c c common /limit/xmin,xmax,ymin,ymax c common /zzhv/hv,plonsw c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine aaxis(xtic,ytic) real xpt,ypt,xmin,xmax,ymin,ymax,xscal,yscal byte hv,plonsw common /limit/xmin,xmax,ymin,ymax common /zzhv/hv,plonsw c********************************************************************* c do some preliminary calculations c********************************************************************* xscal=xmax-xmin yscal=ymax-ymin xpt=xmin !only for vt125 ypt=ymin if (hv) goto 10 xpt=xmin-xscal*0.014 !horizontal axis offset ypt=ymin-yscal*0.014 !vertical axis offset call settl(0.0,1.0) !set tick length c********************************************************************* c now draw axis c********************************************************************* 10 call jump(xmin,ypt) !x axis call draw(xmax,ypt) if(xtic.lt.0.001)goto 40 20 do 30 i=ifix(xmin/xtic),ifix(xmax/xtic) !draw in tics xt=i*xtic call jump(xt,ypt) call ttick(.true.) 30 continue 40 call jump(xpt,ymin) !y axis call draw(xpt,ymax) if(ytic.lt.0.001)goto 70 50 do 60 i=ifix(ymin/ytic),ifix(ymax/ytic) !draw in tics yt=i*ytic call jump(xpt,yt) call ttick(.false.) 60 continue 70 call jump(xpt,ypt) !jump to "zero" call settl(0.5,0.5) !reset default tic length return end c--------------------------------------------------------------------- c c NAME c autosc.for c c USAGE c call autosc(max,min,tic) c c ABSTRACT c Subroutine for automatic scaling of maximum, minimum c and tic limits for axis plotting. c c INPUT c max - maximum vector limit c min - minimum vector limit c c OUTPUT c max - maximum scale limit c min - minimum scale limit c tic - tic interval c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine autosc(min,max,tic) integer iticno real max,min,range real atic,tic,ticno c********************************************************************* c finding the initial tic value c********************************************************************* range=max-min atic=alog10(range) !get order of magnitude if (atic.ge.0) goto 10 !range > 1 atic=atic-1 10 iticno=atic tic=0.5*10**iticno !initial tic size c********************************************************************* c determining the maximum limit c********************************************************************* ticno=max/tic iticno=ticno if (ticno.eq.iticno) goto 30 20 iticno=iticno+1 30 do 40 i=7,13,2 !to get rid of "ugly" no. 40 if (iticno.eq.i) goto 20 !we assume max is +ve if (iticno.le.17) goto 50 iticno=20 50 max=iticno*tic c********************************************************************* c determining the minimum limit c********************************************************************* ticno=min/tic iticno=ticno if (iticno.le.ticno) goto 70 60 iticno=iticno-1 70 k=iabs(iticno) !min might be -ve do 80 i=7,13,2 !to get rid of "ugly" no. 80 if (k.eq.i) goto 60 if (k.le.17) goto 90 k=20 90 iticno=isign(k,iticno) min=iticno*tic c********************************************************************* c determining the tic interval c********************************************************************* 100 range=max-min atic=tic 110 ticno=range/atic if (ticno.gt.3) goto 120 atic=atic/2 !too few tics goto 110 120 if (ticno.lt.9) goto 130 atic=atic+tic !too many tics goto 110 130 tic=atic c********************************************************************* c re-check minimum value c********************************************************************* if (min.le.0) goto 140 !-ve, no worries atic=min-tic if (atic.gt.0) goto 140 !not close to zero min=0.0 goto 100 140 continue c********************************************************************* c the end c********************************************************************* return end c--------------------------------------------------------------------- c c NAME c errbar.for c c USAGE c call errbar(x,y,xerr,yerr) c c ABSTRACT c This subroutine draws in error bars around the data c point. c c INPUT c x,y - point co-ordinates c xerr - error in abscissa vector c yerr - error in ordinate vector c c common /limit/xmin,xmax,ymin,ymax c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine errbar(x,y,xerr,yerr) real x,xerr,xmin,xmax,xpt,xscal real y,yerr,ymin,ymax,ypt,yscal common /limit/xmin,xmax,ymin,ymax c********************************************************************* c some basic calculations c********************************************************************* xscal=(xmax-xmin)*0.015 !unit is 1.5% x range yscal=(ymax-ymin)*0.015 !unit is 1.5% y range c********************************************************************* c draw abscissa error c********************************************************************* if (xscal.ge.xerr) goto 20 !error less than unit xpt=x-xerr if (xpt.lt.xmin) goto 10 !error outside axis? call jump(xpt,y) !goto lower error limit call ttick(.true.) call jump(xpt,y) xpt=x-xscal call draw(xpt,y) !draw lower error bar 10 xpt=x+xscal call jump(xpt,y) xpt=x+xerr if (xpt.gt.xmax) goto 20 !error outside axis? call draw(xpt,y) !draw upper error bar call ttick(.true.) c********************************************************************* c draw ordinate error c********************************************************************* 20 continue if (yscal.ge.yerr) goto 40 !error less than unit ypt=y-yerr if (ypt.lt.ymin) goto 30 !error outside axis? call jump(x,ypt) !goto error limit call ttick(.false.) call jump(x,ypt) ypt=y-yscal call draw(x,ypt) !draw lower error bar 30 ypt=y+yscal call jump(x,ypt) ypt=y+yerr if (ypt.gt.ymax) goto 40 !error outside axis? call draw(x,ypt) !draw upper error bar call ttick(.false.) c********************************************************************* c return to point and program c********************************************************************* 40 continue call jump(x,y) return end c--------------------------------------------------------------------- c c NAME c filnam.for c c ABSTRACT c This subroutine requests a file name from the user c and checks for its' validity. The file may be new, where c the routine checks for any files already existing, or old, c where the routine checks for its' existence. c c INPUT c typ - byte value indicating new (0,f) or old c (1,t) file c input - ascii string, 15 characters long, c format being [dev:filnam.ext]. Arguments c are optional. c ext - ascii string, 4 characters long, default c extension. c c OUTPUT c input - ascii string with file name c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine filnam(typ,input,ext) byte str(69),input(15),ext(4) logical*1 typ c********************************************************************* c get the file name from the user c********************************************************************* 10 type 20 20 format ('$Please type in file name [dev:filnam.ext] > ') accept 30,(input(i),i=1,14) 30 format (14a1) call trim(input) call index(input,':',,m) if (m.ne.0) goto 40 !device specified call concat('dk:',input,str,14) !push in default call scopy(str,input,14) ! device 40 continue call index(input,'.',,m) if (m.ne.0) goto 50 !extension specified call concat(input,'.',input,11) !push in default call concat(input,ext,input,14) ! device 50 continue call scopy('?FILNAM-I-',str,10) !put name in msg str str(11)=27 !esc call concat(str,'[1m',str,14) !bold if (typ) goto 70 c********************************************************************* c check to see if new file already exists c********************************************************************* open(unit=1,type='OLD',name=input,readonly,err=90) close(unit=1,dispose='SAVE') !only if 'new' file call concat(str,input,str) ! is found call concat(str,27,str) call concat(str,'[m already exists, try new file.',str) type 60,7,str 60 format (' ',70a1) goto 10 c********************************************************************* c check to see if existing file is there c********************************************************************* 70 continue open(unit=1,type='OLD',name=input,readonly,err=80) close(unit=1,dispose='SAVE') goto 90 80 continue !only if 'old' file call concat(str,input,str) ! is not found call concat(str,27,str) call concat(str,'[m not on device, try again.',str) type 60,7,str !prompt for new name goto 10 c********************************************************************* c let's go home c********************************************************************* 90 continue return end c--------------------------------------------------------------------- c c NAME c ireply.for c c USAGE c byte=ireply() c if (ireply()) ... c c ABSTRACT c This function gets the input to a request and c determines whether the answer was .true. (y) or .false. c (n). No other replies are allowed. c c INPUTS c ans - byte character c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- logical*1 function ireply() byte ans 10 accept 20,ans 20 format (a1) if (ans.eq.'y') ans='Y' if (ans.eq.'n') ans='N' if (ans.eq.'Y' .or. ans.eq.'N') goto 40 type 30,7 30 format ('$',a1,'Please answer Y or N') goto 10 40 ireply=(ans.eq.'Y') return end c--------------------------------------------------------------------- c c NAME c c labax.for c c USAGE c call labax(idir,rmin,rmax) c c ABSTRACT c Subroutine to label axis tics for GRAPH. It centres c the x axis tic labels under the tic. c c INPUT c idir - direction of axis, 0 = x, 1 = y c rmin - minimum axis value c rmax - maximum axis value c tic - tic interval c c common /limit/xmin,xmax,ymin,ymax c common /zzhv/hv,plonsw c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine labax(idir,rmin,rmax,tic) real rmin,rmax,tic real x,xmin,xmax,xscal,y,ymin,ymax,yscal,f byte hv,plonsw,str(82) integer idir common /limit/xmin,xmax,ymin,ymax common /zzhv/hv,plonsw c********************************************************************* c do some preliminary calculations c********************************************************************* xscal=(xmax-xmin)*0.0431 !1cm x axis (for A4 yscal=(ymax-ymin)*0.0611 !1cm y axis paper) c********************************************************************* c now find length of tic label c********************************************************************* do 80 i=ifix(rmin/tic),ifix(rmax/tic) encode(10,5,str) float(i*tic) 5 format(f10.3) do 10 j=10,7,-1 if (str(j).le."60) str(j)=0 !remove trailing '0''s if (str(j).ne.0) goto 20 !no more 0's, so out 10 continue 20 continue if (idir.eq.0) goto 30 !horizontal axis c********************************************************************* c position y axis tics c********************************************************************* k=len(str) f=0.7 !posn factor for axis if (hv) f=1.0 ! offset x=xmin-xscal*(k*0.2+f) !0.2cm is char width y=i*tic if (hv) goto 70 if (y.eq.ymax) y=ymax-yscal*0.25 !stay within window goto 70 c********************************************************************* c position x axis tics c********************************************************************* 30 continue do 50 l=1,10 if (str(1).gt."40) goto 60 !look for non-sp 1st chr do 40 m=1,10 !not yet so str(m)=str(m+1) ! shift string left 1 chr 40 continue str(10)=0 50 continue 60 continue !now no blank preceeding k=len(str) ! string f=0.2 !chr spacing factor if (hv) f=0.3 ! (ie the width) x=i*tic-k*0.5*xscal*f y=ymin-yscal*0.8 if (hv) y=ymin-yscal*0.4 c********************************************************************* c label axis tics c********************************************************************* 70 continue call labexy(x,y,0.14,0.,str) !char are 0.14cm wide 80 continue ! on the plotter return end c--------------------------------------------------------------------- c c NAME c llabel.for c c USAGE c call llabel(dir,str) c c ABSTRACT c Subroutine to label axis for GRAPH. It writes labels c as appropriate to the axis and centres the text. c c INPUT c dir - label direction (0-horizontal, 1-vertical) c str - string containing label (up to 80 chrs) c c common /limit/xmin,xmax,ymin,ymax c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine llabel(dir,str) integer dir,scal real x,y,xmin,xmax,ymin,ymax byte str(82),hv,plonsw common /limit/xmin,xmax,ymin,ymax common /zzhv/hv,plonsw c********************************************************************* c find out where we want to write the label c********************************************************************* xscal=xmax-xmin !axis ranges yscal=ymax-ymin i=len(str) r=float(i)/2 !centre of string if (dir.eq.0) goto 10 c********************************************************************* c writing veritcal text c********************************************************************* if (hv) x=xmin-xscal*0.1 if (.not.hv) x=xmin-xscal*0.087 y=ymin+yscal*(0.5-r*0.0183) !factor is 3mm y axis call labexy(x,y,0.2,90.,str) goto 20 c********************************************************************* c writing horizontal text c********************************************************************* 10 continue if (hv) y=ymin-yscal*0.065 if (.not.hv) y=ymin-yscal*0.1 x=xmin+xscal*(0.5-r*0.0129) !factor is 3mm x axis call labexy(x,y,0.2,0.,str) 20 return end c--------------------------------------------------------------------- c c NAME c ppoint.for c c USAGE c call ppoint(x,y,ntype) c c ABSTRACT c This subroutine draws a symbol at the present position c depending on the type of point required. c c INPUT c x,y - co-ordinates for point c ntype - integer representing point style c 0 = point 6 = point with circle c 1 = cross 7 = star c 2 = open circle 8 = closed circle c 3 = open triangle 9 = closed triangle c 4 = open square 10= closed square c 5 = open diamond 11= closed diamond c c common /zzhv/hv,plonsw,flh,flv c common /limit/xmin,xmax,ymin,ymax c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine ppoint(x,y,ntype) integer ntype real x,y,xscal,yscal,xpt,ypt,xvt,yvt real xmin,xmax,ymin,ymax byte esc byte hv,plonsw,flh,flv common /zzhv/hv,plonsw,flh,flv common /limit/xmin,xmax,ymin,ymax data esc/27/ c********************************************************************* c work out size and then go to symbol c********************************************************************* xscal=(xmax-xmin)*0.006 !unit is 0.6% x axis yscal=(ymax-ymin)*0.006 !unit is 0.6% y axis ntype=mod(ntype,12) !only 12 sym available call jump(x,y) if (ntype.eq.0 .or. ntype.eq.6) goto 10 !point if (ntype.eq.1 .or. ntype.eq.7) goto 100 !cross if (ntype.eq.2 .or. ntype.eq.8) goto 200 !circle if (ntype.eq.3 .or. ntype.eq.9) goto 300 !triangle if (ntype.eq.4 .or. ntype.eq.10) goto 400 !square if (ntype.eq.5 .or. ntype.eq.11) goto 500 !diamond goto 600 c********************************************************************* c plot a point (without or with a circle) c********************************************************************* 10 continue call draw(x,y) !dot if (ntype.eq.0) goto 600 goto 200 !draws a circle c********************************************************************* c draw a cross (or a star) c********************************************************************* 100 continue !cross xpt=x+xscal ypt=y+yscal call jump(xpt,ypt) !top right corner xpt=x-xscal ypt=y-yscal call draw(xpt,ypt) !bottom left corner xpt=x+xscal call jump(xpt,ypt) !bottom right corner xpt=x-xscal ypt=y+yscal call draw(xpt,ypt) !top left corner if (ntype.eq.1) goto 600 xpt=x+xscal*1.2 !star ypt=y call jump(xpt,ypt) !north xpt=x-xscal*1.2 call draw(xpt,ypt) !south xpt=x ypt=y+yscal*1.2 call jump(xpt,ypt) !east ypt=y-yscal*1.2 call draw(xpt,ypt) !west goto 600 c********************************************************************* c draw a circle (open or closed) c********************************************************************* 200 continue if (ntype.ne.8) goto 230 if (flh) write(9,210)xscal*1.12 !closed circle if (flv) write(7,220)esc 210 format ('WG ',f10.3,',0,360;') 220 format ('+',a1,'Ppw(s1[,])') 230 continue !open circle if (flh) write(9,240)xscal*1.12 !on plotter 240 format ('CI ',f10.3,',5;') if (flv) goto 250 goto 600 250 xpt=x+xscal*1.12 !open circle ypt=y !on screen call scl(xpt,ypt,xvt,yvt) write (7,260)esc,xvt,yvt,esc 260 format ('+',a1,'Ppc(a360)[',f5.0,',',f5.0,']w(s0)',a1,'\') goto 600 c********************************************************************* c draw a triangle (open or closed) c********************************************************************* 300 continue xpt=x-xscal*1.25 !1.25 is optical ypt=y-yscal*0.577*1.25 ! correction factor call jump(xpt,ypt) ! The other no. is if (ntype.ne.9) goto 320 ! for centreing if (flh) write(9,310)xscal*2.5 !closed triangle if (flv) write(7,220)esc 310 format ('WG ',f10.3,',-300,-60,120;') 320 continue !open triangle if (flh) write(9,330)xscal*2.5 330 format ('EW ',f10.3,',-300,-60,120;') if (flv) goto 340 goto 600 340 xvt=xpt+xscal*2.5 yvt=ypt call draw(xvt,yvt) xvt=xpt+xscal*1.25 yvt=ypt+yscal*1.732*1.25 call draw(xvt,yvt) call draw(xpt,ypt) goto 600 c********************************************************************* c draw a square (open or closed) c********************************************************************* 400 continue xpt=x-xscal ypt=y-yscal xvt=x+xscal yvt=y+yscal call jump(xpt,ypt) if (ntype.ne.10) goto 420 if (flh) write(9,410)xvt,yvt !closed square 410 format ('RA ',2f10.3,';') if (flv) write(7,220)esc 420 call draw(xvt,ypt) !open square call draw(xvt,yvt) call draw(xpt,yvt) call draw(xpt,ypt) goto 600 c********************************************************************* c draw a diamond (open or closed) c********************************************************************* 500 continue xpt=x+xscal*1.12 ypt=y+yscal*1.40 xvt=x-xscal*1.12 yvt=y-yscal*1.40 if (ntype.ne.11) goto 540 if (flv) goto 530 !closed diamond call jump(x,yvt) write (9,510)xscal*1.80 510 format ('WG ',f10.3,',51,78,78;') call jump(x,ypt) write (9,520)xscal*1.80 520 format ('WG ',f10.3,',231,78,78;') goto 540 530 write (7,220)esc 540 call jump(x,yvt) !open diamond call draw(xvt,y) call draw(x,ypt) call draw(xpt,y) call draw(x,yvt) goto 600 c********************************************************************* c end of routine c********************************************************************* 600 continue if (ntype.lt.8) goto 620 if (flv) write(7,610)esc,esc !shading off 610 format ('+',a1,'Ppw(s0)',a1,'\') 620 call jump(x,y) return end c--------------------------------------------------------------------- c c NAME c ttick c c USAGE c call ttick(flag) c c ABSTRACT c This subroutine draws a flag at the current position. c A horizontal tick is drawn (ie vertical axis) when the flag is c false and a vertical tick is drawn (ie horizontal axis) when c the flag is true. c c INPUT c flag - logical parameter, T=x axis, F=y axis c c AUTHOR c SYP Wong 28 Jul 86 c c--------------------------------------------------------------------- subroutine ttick(flag) common /zzhv/hv,plonsw byte hv,plonsw common /zzhpb/etx,si,so,esc,EPp(3),Ebsl(2),dum byte etx,si,so,esc,EPp,Ebsl,dum byte flag,char1,char2 c********************************************************************* c check which device c********************************************************************* if (hv) goto 10 if (.not.hv) goto 30 goto 60 c********************************************************************* c write to vt125 c********************************************************************* 10 continue char1='4' !assume horizontal tick (Y axis) char2='0' if(flag)char1='6' !vertical tick (X axis) if(flag)char2='2' type 20,EPp,(char1,i=1,4),(char2,i=1,9),Ebsl 20 format('+',3a1,'P'4a1,'V',9a1,2a1) goto 60 c********************************************************************* c write to hp7475 c********************************************************************* 30 continue if(flag)write(9,40) if(.not.flag)write(9,50) 40 format('XT;') !veritcal X-tick 50 format('YT;') !horizontal Y-tick c********************************************************************* c return c********************************************************************* 60 continue return end