integer function cray2ieg(i) cray2ieg = 1 end c set the pan instruction up subroutine setpan(inst, prange, itype) include 'typecom' c itype = 0 = right, itype = 1 = left, c itype = 2 = up, itype = 3 = down character*(*) inst double precision prange double precision xrange, yrange common/limits/xminp,xmaxp,yminp,ymaxp xrange = (xmaxp - xminp)*prange/200.0d0 yrange = (ymaxp - yminp)*prange/200.0d0 if(itype .eq. 0)then write(inst, *)'dflimits (',xminp+xrange, + yminp, xmaxp+xrange, ymaxp, + '); n; g' else if(itype .eq. 1)then write(inst, *)'dflimits (',xminp-xrange, + yminp, xmaxp-xrange, ymaxp, + '); n; g' else if(itype .eq. 2)then write(inst, *)'dflimits (',xminp, + yminp+yrange, xmaxp, ymaxp+yrange, + '); n; g' else if(itype .eq. 3)then write(inst, *)'dflimits (',xminp, + yminp-yrange, xmaxp, ymaxp-yrange, + '); n; g' else call xtext('Internal error in S/R setpan') end if end c reset loop counter to 1 c (This counter is for the no. of instructions/data lines c stored within each loop) subroutine bloop() include 'typecom' include 'robcom' lnow = 1 end c for storing instructions during a loop subroutine ploop(buff) include 'typecom' character*(*) buff include 'robcom' if(lnow .le. ilmax)then lbuff(lnow) = buff else call xtext('Loop has too much in it') end if lnow = lnow+1 end c for getting stored instructions during a loop subroutine gloop(buff) include 'typecom' character*(*) buff include 'robcom' if(lnow .le. ilmax)then buff = lbuff(lnow) else c no error message - ploop probably does more than enough continue end if lnow = lnow+1 end c decode an instruction into multiple instructions c separated by semi-colons subroutine noinst(inst, instn, ninst) include 'typecom' character*(*) inst character*(*) instn(*) integer ninst include 'charcom' c parameter(iscoln = ichar(';')) integer i, n, ic, istart ninst = 0 n = nnl(inst) istart = 1 c testing 1,2 3 c ninst = 1 c instn(1) = inst(istart:n) c return do 10 i = 1, n ic = ichar(inst(i:i)) if(ic .eq. iscoln)then ninst = ninst + 1 write(instn(ninst), *)inst(istart:i-1) istart = i+1 end if 10 continue ninst = ninst + 1 instn(ninst) = inst(istart:n) end c indentify the nearest point to the specified value c (providing its inside the plotting area) subroutine nneigh(xpos, ypos, xminp,yminp,xmaxp,ymaxp, + x, y, delx, dely) include 'typecom' include 'msizcom' include 'robcom' double precision x(*), y(*), delx(*), dely(*) c only valid if we have data if(npts.le.0)return c only valid in plot area if(xpos.lt.xminp.or.xpos.gt.xmaxp.or. + ypos.lt.yminp.or.ypos.gt.ymaxp) return xscale = (xmaxp - xminp) * (xmaxp - xminp) yscale = (ymaxp - yminp) * (ymaxp - yminp) distmin = (xpos-x(1))*(xpos-x(1))/xscale + + (ypos-y(1))*(ypos-y(1))/yscale imin = 1 do 10 i = 1, npts dist = (xpos-x(i))*(xpos-x(i))/xscale + + (ypos-y(i))*(ypos-y(i))/yscale if(dist .lt. distmin)then distmin = dist imin = i end if 10 continue write(tstring,*)'Nearest point is number',imin call xtext(tstring) write(tstring,*)'x = ',x(imin), ', y = ', y(imin) call xtext(tstring) if(delx(imin) .ne. 0.0d0)then write(tstring,*)'delta x = ', delx(imin) call xtext(tstring) end if if(dely(imin) .ne. 0.0d0)then write(tstring,*)'delta y = ', dely(imin) call xtext(tstring) end if end c set line style to the appropriate type subroutine slstyle(lstyle) include 'typecom' character*(*) lstyle if(lstyle.eq.'SOLID')then call nodash() else if(lstyle.eq.'DASH')then call dodash() else if(lstyle.eq.'DOTDASH')then call dotdash() else if(lstyle.eq.'DOT')then call dodot() else call xtext('ERROR - sltstyle called with unknow style') call xtext(lstyle) end if end c routine to draw either a polyline or a polygon c (difference is whether first point is connected to last or not) c itype specifies which is done subroutine pline(itype, kboard) include 'typecom' integer itype logical kboard include 'robcom' common/inrob/inter logical inter integer press if(.not. kboard)then call xtext('Use cursor to select points') call xtext('Left button selects, middle button terminates') end if call sprompt('Specify start of line') call sprompt('( = use cursor)') ainbuf = ' ' if(kboard.or.(.not.inter).or.HAVEDATA)call getit(ainbuf, 0, icancel) if(icancel .lt. 0)return if(ainbuf.ne.' ')then call savdata(ainbuf) call dcode(ainbuf,ain1,ainfix,aingrd,k) else call cursor(ain1(1), ain1(2)) call fvout(ain1, 2) end if xfirst = ain1(1) yfirst = ain1(2) xtemp = xfirst ytemp = yfirst call movxy(xfirst, yfirst) 102 continue if(kboard.or.(.not.inter).or.HAVEDATA)then call sprompt('Next coordinate, terminates') call getit(ainbuf, 0, icancel) if(icancel .lt. 0)return call dcode(ainbuf,ain1,ainfix,aingrd,k) if(ainbuf.eq.' ')then press = 2 else press = 1 end if else call rband(xtemp, ytemp, 0, ain1(1), ain1(2), 1+2, press) end if if(press.eq.1)then c floating point comparison - to try to save writing duplicate numbers c out if(ain1(1).ne.xtemp.and.ain1(2).ne.ytemp)then call linxy(ain1(1), ain1(2)) xtemp = ain1(1) ytemp = ain1(2) call fvout(ain1, 2) end if goto 102 end if c terminate stored data with a blank line call savdata(' ') c close if polygon if(itype.eq.2)then call linxy(xfirst, yfirst) end if end c draw a filled polygon subroutine pfill(kboard) include 'typecom' logical kboard include 'robcom' common/inrob/inter logical inter integer press call xtext('Use cursor to select points') call xtext('Left button selects, middle button terminates') call sprompt('Specify start of line') call sprompt('( = use cursor)') ainbuf = ' ' if(kboard.or.(.not.inter).or.HAVEDATA)call getit(ainbuf, 0, icancel) if(icancel .lt. 0)return if(ainbuf.ne.' ')then call savdata(ainbuf) call dcode(ainbuf,ain1,ainfix,aingrd,k) else call cursor(ain1(1), ain1(2)) call fvout(ain1, 2) end if xfirst = ain1(1) yfirst = ain1(2) xtemp = xfirst ytemp = yfirst call arkpoly1(xfirst, yfirst) 102 continue if(kboard.or.(.not.inter).or.HAVEDATA)then call sprompt('Next coordinate, terminates') call getit(ainbuf, 0, icancel) if(icancel .lt. 0)return call dcode(ainbuf,ain1,ainfix,aingrd,k) if(ainbuf.eq.' ')then press = 2 else press = 1 end if else call rband(xtemp, ytemp, 0, ain1(1), ain1(2), 1+2, press) end if if(press.eq.1)then c floating point comparison - to try to save writing duplicate numbers c out if(ain1(1).ne.xtemp.and.ain1(2).ne.ytemp)then c we indicate where the polygon is going to be using lines c however, switch off storing these in the "postscript" array c as this isn't needed call haltps() call linxy(ain1(1), ain1(2)) call startps() call arkpoly(ain1(1), ain1(2)) xtemp = ain1(1) ytemp = ain1(2) call fvout(ain1, 2) end if goto 102 end if c terminate stored data with a blank line call savdata(' ') c terminate polygon call arkpoly2() end c call the ruler routine subroutine ruler() include 'typecom' include 'robcom' double precision x1, y1, x2, y2, r integer press 10 continue c call cursor(x1, y1) call bcursor(x1, y1) call rulerp(x1, y1, x2, y2, 5, press) if(press .eq. 3)goto 10 write(tstring, *)'Coordinates were: (',x1,',',y1, '), (', + x2,',',y2,')' call xtext(tstring) write(tstring, *)'Center value (x, y) = (', (x1+x2)/2, ',', + (y1+y2)/2, ')' call xtext(tstring) write(tstring, *)'Delta X =', x2-x1, ', Delta Y = ', y2-y1 call xtext(tstring) r = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1)) write(tstring, *)'Total distance = ', r call xtext(tstring) end c general purpose call to rubber banding routines subroutine nline(x1, y1, x2, y2, itype, s1, s2) include 'typecom' integer press character*(*) s1, s2 press = 0 call update() 10 call xtext(s1) call cursor(x1, y1) call xtext(s2) call rband(x1, y1, itype, x2, y2, 5, press) if(press .eq. 3) goto 10 end c getnum - a call to getit followed by a call to decode c part of a slow change to improved handling of input subroutine getnum(doit, ain, ifix, igrd, k) include 'typecom' integer doit double precision ain(*) integer ifix(*), igrd(*) integer k c a single common block from the robcom file include 'robcom' call getit(ainbuf, doit, icancel) if(icancel .lt. 0)then k = icancel return end if call dcode(ainbuf, ain, ifix, igrd, k) end c requests a specified number of arguments subroutine getinum(doit, ain, ifix, igrd, k, iwant) include 'typecom' common/inrob/inter logical inter integer doit double precision ain(*) integer ifix(*), igrd(*) integer iwant integer k c a single common block from the robcom file include 'robcom' 10 call getit(ainbuf, doit, icancel) c test for if the "cancel" button was pressed if(icancel .lt. 0)then k = icancel return end if call dcode(ainbuf, ain, ifix, igrd, k) if(k .ne. iwant)then call sprompt('Wrong number of parameters') write(tstring, *)iwant,' needed,', k, ' given' call sprompt(tstring) if(inter)then goto 10 else call xtext('WARNING: Wrong number of parameters read') end if end if end c getit is intended to be the main I/O routine for c dialogue interaction within Robot. c It will either read from the screen or from a file c depending on whether INTER is true or not. c This now also automatically writes the input to unit c no. 7 too if DOIT is 1. subroutine getit(buff, doit, iscancel) include 'typecom' common/inrob/inter common/stdirob/stdin logical stdin logical inter include 'robcom' character*(*) buff character*200 bsave integer doit integer iscancel logical cancel c obtain a text string, either read from a file c if non-intereactive of else call up the appropriate X c windows subroutine c first force the cancel option to be inoperative iscancel = 1 c force increase of getittime to try and keep menus in step... call incgt() c first option is to check whether we already have data from c somewhere else DONEDATA = .FALSE. bsave = buff buff = ' ' if(HAVEDATA)then DONEDATA = .TRUE. read(bufinst(ibufstr),'(a)', err = 900)buff ibufstr = ibufstr + 1 c see whether we've now run out of stored data if(ibufstr .ge. maxbuf)HAVEDATA = .FALSE. c else if(stdin)then else if(stdin .and. inter)then read(*, '(a)',err =900,end=901)buff else if(inter) then iscancel = 1 call gettextxv(buff, iscancel) c restore string if cancelled if(cancel(iscancel))buff = bsave call blankit(buff) else if((.not.(loop)).or.rloop.eq.1) then read(ifil, '(a)', err = 900, end = 900)buff else if(loop.and.rloop.gt.1)then if(ifil .eq. uloop)then call gloop(buff) else read(ifil, '(a)',err =900,end=901)buff end if end if end if c C routines leave a NULL at the end of a character string c to make sure there is no junk after this change everything c after a null to a space call blankit(buff) c how data is stored depends on value of doit and whether we're c looping. If looping just store data for first pass for retrieval c on subsequent passes if(doit.eq.1) call savdata(buff) c reset prompt call preset() return 900 call xtext('Error occurred during a read from the input file') return 901 call xtext('Internal Error: End of file detected in getit') end c save data in unit number 7 and/or write in arrays used in looping c as required subroutine savdata(buffin) include 'typecom' c character buffin(*) character*(*) buffin character*200 buff include 'robcom' integer imax, imaxin, itodo, i imaxin = len(buffin) imax = len(buff) c print*,'imax=', imax, 'imaxin = ',imaxin c clean buff do 5 i = 1, imax buff(i:i) = ' ' 5 continue c set as much of buff equal to buffin as possible itodo = MIN(imax, imaxin) do 10 i = 1, itodo buff(i:i) = buffin(i:i) 10 continue call blankit(buff) if(.not.(loop))then c write(ifil7,'(a)')buff(1:nnl(buff)) call pseudo7(buff, 1) end if if(loop.and.rloop.eq.1)then c write(ifil7,'(a)')buff(1:nnl(buff)) call pseudo7(buff, 1) if(.not.(DONEDATA) .and. + ifil .eq. uloop)then call ploop(buff) end if end if DONEDATA = .FALSE. end c pseudo7 is a replacement for writing directly to file 7 c idea is to store output for file 7. This is to try to enable c "cancel" options to work without leaving broken stuff in c log file. c This subroutine either stores an instruction in a (limited) buffer, c flushes it out, or writes it to the real file 7. Controlled by c value of "op" input parameter subroutine pseudo7(textin, op) include 'typecom' character*(*) textin integer op save include 'robcom' c local variables integer lmax, linmax parameter(lmax = 15) parameter(linmax = 200) character*(linmax) blocal(lmax) integer lineno integer i c reinitalize local storage for op == 0 if(op .eq. 0)then lineno = 1 c add an instruction for op == 1 else if(op .eq. 1)then if(lineno .le. lmax)then blocal(lineno) = textin lineno = lineno + 1 else do 5 i = 1, lmax tstring = blocal(i) write(ifil7,'(a)')tstring(1: nnl(tstring)) 5 continue c and write out the lastest unbuffered instruction write(ifil7,'(a)')textin(1: nnl(textin)) c and reset line indicator lineno = 1 c print*,'ERROR IN PSUEDO7, BUFFER TOO SMALL' c print*, 'lineno = ', lineno c print*,'Data was' c print*,textin(1:nnl(textin)) end if c print existing instructions for op == 2 else if(op .eq. 2)then if(lineno .eq. 1)return do 20 i = 1, MIN((lineno - 1), lmax) tstring = blocal(i) write(ifil7,'(a)')tstring(1: nnl(tstring)) 20 continue c and reset line indicator lineno = 1 c catch unrecognized op value and print error message else print*,'ERROR - UNEXPECTED INPUT TO PSEUDO7' print*,'op = ', op end if end c in a kludgey type of way see whether the instruction has anything c after the main instruction - if so assume it's data subroutine chkdat(inst) include 'typecom' character*(*) inst character*200 buffer imax = len(inst) buffer = ' ' c find where first blank occurs do 10 i = 1, imax if(inst(i:i).eq.' ')then goto 11 end if 10 continue c didn't find any blank I think return 11 continue j = i c check whether there's a NON-blank after this do 20 i = j, imax if(inst(i:i).ne.' ')then goto 21 end if 20 continue c didn't find any non-blank I suppose return 21 continue buffer = inst(i:imax) call buffset(buffer) c and finally remove the nasty old data from the instruction do 30 i = j, imax inst(i:i) = ' ' 30 continue end c list the aliases available subroutine lalias(useral1, useral2, noual) include 'typecom' character*(*) useral1(*), useral2(*) include 'robcom' if(noual.eq.0)then call xtext('No aliases in use') else call xtext('Defined aliases are:') end if do 10 i = 1, noual tstring = useral1(i) call xtext(tstring) tstring = useral2(i) call xtext(tstring) 10 continue end c add another alias to the list or change an existing definition subroutine adalias(buff1, buff2, useral1, useral2, noual, maxual) include 'typecom' character*(*) buff1, buff2, useral1(*), useral2(*) c convert to upper case call lctouc(buff1) c now see whether the alias is already in our list do 10 i = 1, noual if(buff1(1:nnl(buff1)).eq.useral1(i))then useral2(i) = buff2 goto 20 end if 10 continue c seems we need to add a new alias to the list if(noual.ge.maxual)then call xtext('ERROR: TOO MANY ALIASES') else noual = noual + 1 useral1(noual) = buff1 useral2(noual) = buff2 end if 20 continue end c check if an instruction is an alias and remove it if it is subroutine chkal(inst, useral1, useral2, noual) include 'typecom' character*(*) inst, useral1(*), useral2(*) character*120 inst2 integer ilen inst2 = inst call lctouc(inst2) ilen = nnl(inst2) do 10 i = 1, noual if(inst2(1:ilen).eq.useral1(i))then inst = useral2(i) c an additional data check in addition to that in the main routine c call chkdat(inst) c and make it upper case c call lctouc(inst) goto 20 end if 10 continue 20 continue end c remove a name from the alias list subroutine dalias(buff, useral1, useral2, noual) include 'typecom' character*(*) buff, useral1(*), useral2(*) do 10 i = 1, noual if(buff(1:nnl(buff)).eq.useral1(i))then do 20 j = i, noual useral1(j) = useral1(j+1) useral2(j) = useral2(j+1) 20 continue noual = noual - 1 goto 30 end if 10 continue call xtext('ERROR: ALIAS NOT FOUND') 30 continue end c used by external routines to set up a buffer which can be read c from within robot subroutine buffset(bufft) include 'typecom' character*120 buffer character bufft(*) include 'robcom' integer imax, i imax = len(buffer) do 10 i = 1, imax buffer(i:i) = bufft(i) 10 continue c check if there's enough space in the buffer if(ibufstr.gt.ibufmax)then print*,'ERROR: not enough space in buffers.' print*,'complain to corbet@lheamail.gsfc.nasa.gov' return end if HAVEDATA = .TRUE. call blankit(buffer) write(bufinst(ibufstr),'(a)')buffer ibufstr = ibufstr + 1 maxbuf = ibufstr end c reset data buffer subroutine bufrst() include 'typecom' include 'robcom' HAVEDATA = .FALSE. ibufstr = 1 end c send a value to log file (also potentially stores in loop array via c call to savdata) c fvout - FORTRAN values out subroutine fvout(value, n) include 'typecom' double precision value(*) integer n character*120 mybuff write(mybuff, *)(real(value(i)), i = 1, n) call savdata(mybuff) end c fiout - FORTRAN integer values out subroutine fiout(ivalue, n) include 'typecom' integer ivalue(*) integer n character*120 mybuff write(mybuff, *)(ivalue(i), i = 1, n) call savdata(mybuff) end c send an ***instruction*** file 7 c Hence no call to savdata c faout - FORTRAN ASCII out subroutine saveinst(texti) include 'typecom' character*(*) texti character*120 text include 'robcom' integer imax, i imax = len(text) do 10 i = 1, imax text(i:i) = texti(i:i) 10 continue call blankit(text) if((.not.(loop)).or.rloop.eq.1)then c write(ifil7, '(a)')text(1:nnl(text)) c save in small internal buffer call pseudo7(text, 1) end if c first passage through loop store instruction c if(loop.and.rloop.eq.1)then c write(rinst(rinstno),'(a)')text c rinstno = rinstno + 1 c end if end c This takes a string returned by a C routine and replaces c anything after a NULL with a space c subroutine blankit(buff) c character*(*) buff c integer i, imax c logical isend c isend = .FALSE. c c imax = len(buff) c c do 10 i = 1, imax c if(buff(i:i).eq.char(0)) isend = .TRUE. c if(isend) buff(i:i) = ' ' c10 continue c c end c A faster way? c This takes a string returned by a C routine and replaces c anything after a NULL with a space subroutine blankit(buff) include 'typecom' character*(*) buff integer i, imax imax = len(buff) do 10 i = 1, imax if(buff(i:i) .eq. char(0)) goto 100 10 continue return 100 continue do 110 j = i, imax buff(j:j) = ' ' 110 continue end c test whether we are to carry out the next instruction c because it's an IF related instruction or the c IF tests in operation will let us c returns 0 if OK, otherwise 1 integer function ifcheck(inst, condition, iflevel) character*(*) inst logical condition(*) integer iflevel ifcheck = 0 if(iflevel.eq.0)return if(inst.eq.'IF'.or. + inst.eq.'ELSE'.or. + inst.eq.'ELSEIF'.or. + inst.eq.'ENDIF')return do 10 i = 1, iflevel if(.not.condition(i))then ifcheck = 1 return end if 10 continue end c pause using ca calculation (hence system speed dependent) double precision function waiter(a) n = nint(a) k = n/2 waiter = 0.0d0 do 20 m = 1, k do 10 i = 1, n x = dble(n) y = x**1.2d0 z = cos(x)*sin(y) z2 = sqrt(abs(dble(i-m))) waiter = waiter+z2+z 10 continue 20 continue end c list user variables and their values subroutine lvars() include 'typecom' include 'robcom' do 10 i = 1, nouvar write(ainbuf, *)uservn(i) write(tstring, *)ainbuf(1:nnl(ainbuf)),' =',uservv(i) call xtext(tstring) 10 continue write(tstring, *)'A1 =',auser1 call xtext(tstring) write(tstring, *)'A2 =',auser2 call xtext(tstring) write(tstring, *)'A3 =',auser3 call xtext(tstring) write(tstring, *)'A4 =',auser4 call xtext(tstring) write(tstring, *)'A5 =',auser5 call xtext(tstring) end c add a variable to the existing list subroutine addvar(buff) include 'typecom' character*(*) buff include 'robcom' if(nouvar.ge.maxuvar)then call xtext('ERROR: TOO MANY VARIABLES') write(tstring, *)'MAXIMUM NUMBER = ',maxuvar call xtext(tstring) return end if c it was too complicated to keep case sensitivity call lctouc(buff) c check whether variable already exists do 10 i = 1, nouvar if(uservn(i) .eq. buff)then write(tstring, *)'ERROR: VARIABLE "', + buff(1:nnl(buff)), '" ALREADY IN USE' call xtext(tstring) return end if 10 continue nouvar = nouvar + 1 uservn(nouvar) = buff end c assign a value to a user variable c either a pregrogrammed one (a1 to a5) c or a user specified one subroutine assign(name, value, ifail) include 'typecom' include 'robcom' character*(*) name double precision value integer ifail ifail = 0 do 10 i = 1, nouvar if(name.eq.uservn(i))then uservv(i) = value return end if 10 continue if(name.eq.'A1')then auser1 = value return else if(name.eq.'A2')then auser2 = value return else if(name.eq.'A3')then auser3 = value return else if(name.eq.'A4')then auser4 = value return else if(name.eq.'A5')then auser5 = value return end if ifail = 1 end c simply a stop statement to be called by C! subroutine fstop() include 'typecom' stop end c flush out log and info files subroutine liflsh() include 'typecom' include 'robcom' call flush(ifil7) call flush(ifil8) end c opens a "FILE" containing Robot instructions c if the name passed isn't found we also try adding ".rob" c to the end of it c ifil is set equal to ulevel and then ulevel c is incremented subroutine fopen(newfil, ulevel, ifil, ifail) include 'typecom' character*(*) newfil integer ifail, ulevel, ifil ifail = 0 open(ulevel, file=newfil, status='old', err = 10) c increase ulevel for next read to a read to file ifil = ulevel ulevel = ulevel + 1 return 10 continue newfil = newfil(1:nnl(newfil))//'.rob' open(ulevel, file=newfil, status='old', err = 20) c increase ulevel for next read to a read to file ifil = ulevel ulevel = ulevel + 1 return 20 continue ifail = 1 end integer function nnl(buff) c returns lengths of string which is not spaces or nulls character*(*) buff include 'charcom' c parameter (ispace = ichar(' ')) integer ic ilen = len(buff) do 10 nnl = ilen, 1, -1 ic = ichar(buff(nnl:nnl)) if(ic .ne. ispace .and. ic .ne. 0) goto 11 10 continue 11 continue end c close down units 7 and 8 subroutine clf() include 'typecom' include 'robcom' close(unit = ifil7, err = 1000) close(unit = ifil8, err = 1000) return 1000 call xtext('ERROR: closing log/info file') end c call to create log and info files c called once at start of robot c also called if we delete the log files and want new ones subroutine olf() include 'typecom' include 'robcom' integer ifail ainbuf = 'rl' ifil7 = 7 ifail = 0 call ofile(ifil7, ainbuf, ifail) c store log file name for later use logf = ainbuf c call C routine to set mode to executable and also to c add a line at the top which should enable auto execution if c location of robot was specified correctly c Don't do this if the file couldn't be opened if(ifail .ne. 1) call setfl(ainbuf, n) tstring = '#!'//ainbuf(1:MIN(n, nnl(ainbuf)))//'robot' write(ifil7, '(a)') tstring(1:nnl(tstring)) ainbuf = 'ri' ifil8 = 8 ifail = 0 call ofile(ifil8, ainbuf, ifail) c store info file name for later use infof = ainbuf end c open units 7 (ROBOTLOG) and 8 (ROBOTINF) c if these fail write to /dev/null subroutine ofile(ifile, text, ifail) include 'typecom' include 'robcom' integer ifile, ifail integer fpid character*(*) text character*5 pbuf ipid = fpid() write(pbuf, '(i5.5)')ipid tstring = '/tmp/'//pbuf(1:nnl(pbuf))//text open(ifile, file = tstring, status = 'new', err = 6) text = tstring(1:nnl(tstring)) return 6 print*,'**WARNING***' print*,'COULDN''T OPEN ', tstring print*,'WRITING TO /dev/null' ifail = 1 open(ifile, file = '/dev/null', status = 'old', err = 7) return 7 print*,'COULDN''T EVEN OPEN /dev/null!!!' return end c get log file name - used in C routines subroutine getfl(buff) include 'typecom' include 'robcom' character*(*) buff integer i, length length = nnl(logf) do 10 i = 1, length buff(i:i) = logf(i:i) 10 continue buff(length+1:length+1) = char(0) end c get info file name - used in C routines subroutine getfi(buff) include 'typecom' include 'robcom' character*(*) buff integer i, length length = nnl(infof) do 10 i = 1, length buff(i:i) = infof(i:i) 10 continue buff(length+1:length+1) = char(0) end c write a character string to ifil8 (done this c way so that a call to nnl doesn't have to be written c every time) subroutine write8(string) include 'typecom' character*(*) string include 'robcom' write(ifil8, *)string(1 : nnl(string)) end c read a Ginga spectrum (48 channels) subroutine specred(pfnam,title,xlab,ylab, +x,y,delx,dely,npts) include 'typecom' include 'robcom' double precision x(*),y(*),delx(*),dely(*) character*(*) title,xlab,ylab,pfnam open(unit=42,file=pfnam,status='old',err=903) npts=48 read(42,'(a)')title do 10 i=1,19 10 read(42,'(a)')xlab c get the bin energies read(42,66)itemp,(x(i),i=1,12) read(42,66)itemp,(x(i),i=13,24) read(42,66)itemp,(x(i),i=25,36) read(42,77)itemp,(x(i),i=37,49) read(42,66)itemp,(y(i),i=1,12) read(42,66)itemp,(y(i),i=13,24) read(42,66)itemp,(y(i),i=25,36) read(42,66)itemp,(y(i),i=37,48) read(42,66)itemp,(dely(i),i=1,12) read(42,66)itemp,(dely(i),i=13,24) read(42,66)itemp,(dely(i),i=25,36) read(42,66)itemp,(dely(i),i=37,48) 66 format(i1,12d12.5) 77 format(i1,13d12.5) c sort out x values and errors do 20 i=1,48 delx(i)=(x(i+1)-x(i))/2.0d0 x(i)=(x(i)+x(i+1))/2.0d0 20 continue xlab='Energy (keV)' ylab='Counts' close(unit=42) return 903 call xtext('ERROR IN ROUTINE SPECRED!!!') end subroutine filred(pfnam,title,xlab,ylab,x,y,delx,dely,msiz,npts, +xfit,yfit,xfite,yfite,nptsf,ifail,istar) include 'typecom' include 'robcom' c maxcol is maximum number of columns that filred can cope with parameter (maxcol=25) common/table/index,ftext logical ftext common/echoc/echo logical echo integer index(4) character*(*) title,xlab,ylab,pfnam c arrays for read routine double precision ain(maxcol) integer ifix(maxcol),igrd(maxcol) double precision x(*),y(*),dely(*),delx(*) double precision xfit(*),yfit(*),xfite(*),yfite(*) integer igood ifail=0 igood = 1 open(unit=42,file=pfnam,status='old',err=2) goto 3 2 continue c we failed to open the name specified - let's see if adding c a ".dat" to the end will help us pfnam = pfnam(1:nnl(pfnam))//'.dat' open(unit=42,file=pfnam,status='old',err=903) 3 continue C READ IN DATA FROM PLOT FILE if(ECHO)then write(tstring,*)'Data file name is: ',pfnam call totext(tstring) call xtext('Reading data...') end if c text at top of file? if (ftext)then c identifying text read(42,'(a)',end=994,err=995)title c x title read(42,'(a)',end=994,err=995)xlab c y title read(42,'(a)',end=994,err=995)ylab else c title is file name title = pfnam c and clear x and y labels xlab = ' ' ylab = ' ' end if c zero ain do 4 i=1,maxcol 4 ain(i)=0.0d0 c read in data values i=istar 5 read(42,'(a)',end=50,err=993)ainbuf c Use dcode2 - like dcode but with far fewer if statements c so can't use key words or comments or split variables up c with different types of brackets call dcode2(ainbuf,ain,ifix,igrd,k) c ignore blank lines if(k .eq. 0)goto 5 c if too much data just keep count of how many points if(igood .eq. 0)then i = i + 1 goto 5 end if if(i .gt. msiz)goto 992 if(k.eq.1)then y(i)=ain(1) x(i)=dble(i) delx(i)=0.0d0 dely(i)=0.0d0 else if(index(1) .gt. 0)then x(i) = ain(index(1)) else x(i) = dble(i) end if if(index(2) .gt. 0)then y(i) = ain(index(2)) else y(i) = dble(i) end if if(index(3) .gt. 0)then dely(i) = ain(index(3)) else dely(i) = 0.0d0 end if if(index(4) .gt. 0)then delx(i) = ain(index(4)) else delx(i) = 0.0d0 end if end if xfit(i) = x(i) yfit(i) = y(i) xfite(i) = delx(i) yfite(i) = dely(i) i=i+1 c zero ain do 44 j = 1, k 44 ain(j) = 0.0d0 goto 5 50 continue if(igood .eq. 1)then npts = i - 1 nptsf = npts end if if(echo)then write(tstring,*)'number of points read = ',npts call totext(tstring) end if close(unit=42) c info message if too much data if(igood .eq. 0)then write(tstring, *)'TOTAL SIZE OF INPUT DATA = ',i-1 call totext(tstring) call xtext('USE AT LEAST THAT SIZE WITH "ARRAYSIZE"') end if return C ERROR MESSAGES 992 call xtext('ERROR: THERE IS TOO MUCH DATA') write(tstring,*)'A MAXIMUM OF',MSIZ,' POINTS IS ALLOWED' call totext(tstring) call totext('TO CHANGE THIS USE THE "ARRAYSIZE" COMMAND') call totext('Checking size of input file...') npts = i - 1 nptsf = npts igood = 0 goto 5 903 call xtext('ERROR OPENING DATA FILE') write(tstring, *)'file: ',pfnam(1:nnl(pfnam)-4) call xtext(tstring) IFAIL=1 return 993 call xtext('ERROR READING DATA FROM DATA FILE') write(tstring,*)'AT LINE NUMBER',I call totext(tstring) return 994 call xtext('ERROR: THERE IS NO DATA IN THE DATA FILE') ifail = 1 return 995 call xtext('ERROR READING TEXT HEADER FROM THE DATA FILE') ifail = 1 return end C ENTER DATA DIRECTLY INTO PROGRAM subroutine type(x,y,delx,dely,npts, +title,xlab,ylab,msiz,istar, ilabels) include 'typecom' include 'robcom' C MAXCOL IS MAXIMUM NUMBER OF COLUMNS THAT FILRED CAN COPE WITH parameter (maxcol=25) common/table/index,ftext logical ftext integer index(4) character*(*) title,xlab,ylab C ARRAYS FOR READ ROUTINE double precision ain(maxcol) integer ifix(maxcol),igrd(maxcol) double precision x(*),y(*),dely(*),delx(*) integer igood character*5 helper(4) common/inrob/inter logical inter logical cancel data helper/'x','y','yerr ','xerr '/ igood = 1 C Text at top of file? IF (FTEXT.and.(ilabels.eq.1))then C IDENTIFYING TEXT call sprompt('Give title for graph') call getit(title, 1, icancel) if(icancel .lt. 0)return C X TITLE call sprompt('Give label for X-axis') call getit(xlab, 1, icancel) if(icancel .lt. 0)return C Y TITLE call sprompt('Give label for Y-axis') call getit(ylab, 1, icancel) if(icancel .lt. 0)return end if C ZERO AIN do 4 i=1,maxcol 4 ain(i)=0.0d0 c explain which columns are which write(tstring,*)helper(index(1)),helper(index(2)), + helper(index(3)),helper(index(4)) call sprompt(tstring) C READ IN DATA VALUES i=istar 5 if(inter) call sprompt('Give data, blank line to end') call getnum(1, ain, ifix, igrd, k) if(k .le. 0)goto 50 c flush the pseduo7 buffer - this means that cancel won't really work! c call pseudo7(ainbuf, 2) c Try autoflush in pseudo7 itself c if too much data just keep count of how many points if(igood .eq. 0)then i = i + 1 goto 5 end if if(i.gt.msiz)goto 992 if(k.eq.1)then y(i) = ain(1) x(i) = dble(i) delx(i) = 0.0d0 dely(i) = 0.0d0 else if(index(1) .gt. 0)then x(i) = ain(index(1)) else x(i) = dble(i) end if if(index(2) .gt. 0)then y(i) = ain(index(2)) else y(i) = 0.0d0 end if if(index(3) .gt. 0)then dely(i)=ain(index(3)) else dely(i) = 0.0d0 end if if(index(4) .gt. 0)then delx(i)=ain(index(4)) else delx(i)=0.0d0 end if end if i=i+1 C ZERO AIN do 44 j = 1, k 44 ain(j) = 0.0d0 goto 5 50 continue if(cancel(icancel))then call xtext('WARNING - Cancel limited for this command') call xtext('due to buffer space limitations') c icancel = 0 end if if(igood .eq. 1) npts = i - 1 write(tstring,*)'number of points read = ',npts call totext(tstring) c info message if too much data if(igood .eq. 0)then write(tstring, *)'TOTAL SIZE OF INPUT DATA = ',i-1 call totext(tstring) call xtext('USE AT LEAST THAT SIZE WITH "ARRAYSIZE"') end if return C ERROR MESSAGES 992 write(tstring,'(a, i)')'ERROR: THERE IS TOO MUCH DATA IN UNIT', IFIL call totext(tstring) write(tstring,*)'A MAXIMUM OF',MSIZ,' POINTS IS ALLOWED' call totext(tstring) call totext('TO CHANGE THIS USE THE "ARRAYSIZE" COMMAND') igood = 0 goto 5 c995 write(tstring,'(a, i)')'ERROR READING TEXT HEADER FROM UNIT', IFIL c call totext(tstring) c return end C OUTPUT LABELS AND X,Y POINTS TO SPECIFIED FILE subroutine filwrt(pfnam,title,xlab,ylab,x,y,delx,dely, +npts, clobber) include 'typecom' include 'robcom' common/test/xtmin, xtmax, ytmin, ytmax character*(*) title,xlab,ylab,pfnam double precision x(*), y(*), dely(*), delx(*) integer npts logical clobber if(pfnam.ne.' ')then if(clobber) then open(unit=42, file=pfnam, status = 'unknown', err = 999) else open(unit=42, file=pfnam, status = 'new', err = 999) end if c write data to plot file c identifying text write(42,'(a)')title(1:nnl(title)) c x title write(42,'(a)')xlab(1:nnl(xlab)) c y title write(42,'(a)')ylab(1:nnl(ylab)) c write out data values do 5 i=1,npts if(x(i) .gt. xtmax .or. x(i) .lt. xtmin .or. + y(i) .gt. ytmax .or. y(i) .lt. ytmin) goto 5 write(42,*, err = 1000)x(i), y(i), dely(i), delx(i) 5 continue close(unit=42) else call xtext(title) call xtext(xlab) call xtext(ylab) C WRITE OUT DATA VALUES do 15 i=1,npts if(x(i).gt.xtmax.or.x(i).lt.xtmin.or. + y(i).gt.ytmax.or.y(i).lt.ytmin) goto 15 write(tstring,*)x(i),y(i),dely(i),delx(i) call totext(tstring) 15 continue end if return c error messages 999 call xtext('ERROR: OPENING FILE FOR OUTPUT') call xtext('MAYBE FILE ALREADY EXISTS?') call xtext('- USE "CLOBBER" COMMAND TO ENABLE OVERWRITING') return 1000 call xtext('ERROR: WRITING TO OUTPUT FILE') end C OUTPUT ONLINE SHORT VERSION OF MANUAL DURING INTERACTIVE SESSION subroutine manual include 'typecom' include 'robcom' open(unit=42,file='manual', status = 'old', err=900) 5 read(42,'(a)',end=10, err= 900)ainbuf call xtext(ainbuf) goto 5 10 continue close(unit=42) return 900 call xtext('*** ERROR OPENING MANUAL ***') end subroutine dmax(x,y,dely,delx,xmax,ymax,xmin,ymin,npts) include 'typecom' include 'robcom' c find data limits double precision x(*),y(*),dely(*),delx(*) double precision xdel, ydel c taking absolute values allows for having negative c errors (these may result if the data values are multiplied c by a negative number) xdel = abs(delx(1)) ydel = abs(dely(1)) xmax=x(1)+xdel xmin=x(1)-xdel ymin=y(1)-ydel ymax=y(1)+ydel do 10 i=2,npts xdel = abs(delx(i)) ydel = abs(dely(i)) xmax=max(xmax,x(i)+xdel) ymax=max(ymax,y(i)+ydel) xmin=min(xmin,x(i)-xdel) ymin=min(ymin,y(i)-ydel) 10 continue end c plot data as a filled histogram subroutine fhist(npts, x, y) include 'typecom' double precision x(*), y(*) COMMON/TEST/xtmin, xtmax, ytmin, ytmax C FIND START POINT INSIDE DATA LIMITS DO 18 I=1,NPTS xp = x(i) yp = y(i) if(xp.gt.xtmin.and.xp.lt.xtmax.and.yp.gt.ytmin + .and.yp.lt.ytmax)then c CALL MOVXY(X(I),Y(I)) x1 = x(i) y1 = y(i) GOTO 19 end if 18 continue 19 J=I C EXTEND LINE BACKWARD IF POSSIBLE if(j.gt.1)xp=(x(i)+x(i-1))/2.0d0 if(j.eq.1)xp=x(i)-((x(i+1)-x(i))/2.0d0) yp=y(i) if(xp.lt.xtmin)xp=xtmin if(yp.lt.ytmin)yp=ytmin if(xp.gt.xtmax)xp=xtmax if(yp.gt.ytmax)yp=ytmax call fbox(x1, ytmin, xp, yp) x1 = x(j) y2 = y(j) do 110 i=j+1,npts xp=(x(i)+x(i-1))/2.0d0 yp=y(i-1) if(yp.gt.ytmax)yp=ytmax if(xp.lt.xtmin)xp=xtmin if(yp.lt.ytmin)yp=ytmin if(xp.gt.xtmax)goto 111 call fbox(x1, ytmin, xp, yp) yp=y(i) if(yp.gt.ytmax)yp=ytmax if(yp.lt.ytmin)yp=ytmin x1 = xp xp=x(i) if(xp.gt.xtmax)goto 111 call fbox(x1, ytmin, xp, yp) x1 = xp 110 continue C TERMINATE HISTOGRAM xp=x(npts)+((x(npts)-x(npts-1))/2.0d0) if(xp.gt.xtmax)xp=xtmax call fbox(x1, ytmin, xp, yp) goto 112 111 continue call fbox(x1, ytmin, xtmax, yp) 112 continue return end C PLOT A SYMBOL C The symbol is the ASCII character specified by SYMB c this is not very accurately centred in the Y direction subroutine tplot(npts,x,y,symb) include 'typecom' include 'robcom' common/chsiz/xchsiz, ychsiz, fsize common/test/xtmin, xtmax, ytmin, ytmax double precision x(*),y(*) character*(*) symb logical hasbs hasbs = .FALSE. imax = nnl(symb) c check for control character call getlen(symb, imax, iret, hasbs, xlen) ainbuf = symb do 10 i=1,npts if(x(i).lt.xtmin.or.x(i).gt.xtmax.or.y(i).lt.ytmin.or.y(i).gt. +ytmax)goto 10 if(hasbs)then start=x(i)-xlen/2.0d0 call txtm(start,y(i)-ychsiz/2.3d0,symb,imax) else call style(x(i), y(i)-ychsiz/2.3d0, ainbuf, nnl(ainbuf), 0) end if 10 continue end c plot text centred in x direction on given position subroutine ctext(xpos, ypos, text) include 'typecom' double precision xpos, ypos character*(*) text include 'robcom' common/chsiz/xchsiz, ychsiz, fsize logical hasbs integer imax imax = MIN(len(ainbuf), nnl(text)) hasbs = .FALSE. ainbuf = text(1: imax) c is there a backspace? call getlen(text, imax, iret, hasbs, xlen) if(hasbs)then start=xpos-xlen/2.0d0 call txtm(start,ypos-ychsiz/2.3d0,text,imax) else call style(xpos, ypos-ychsiz/2.3d0, ainbuf, nnl(ainbuf), 0) end if end c plot a symbol c tplot2 plots a limited set of symbols accurately centred c on the point positions (both on screen and in the postscript c outputlines subroutine tplot2(npts,x,y,isymb) include 'typecom' include 'robcom' common/test/xtmin, xtmax, ytmin, ytmax double precision x(*),y(*) integer isymb do 10 i=1,npts if(x(i).lt.xtmin.or.x(i).gt.xtmax.or.y(i).lt.ytmin.or.y(i).gt. +ytmax)goto 10 call symbol(x(i),y(i), isymb) 10 continue end c plot data as pillars (this is more useful in 3D mode) subroutine pillar(npts, x, y) include 'typecom' double precision x(*), y(*) common/test/xtmin, xtmax, ytmin, ytmax do 10 i = 1, npts if(x(i).gt.xtmax.or.x(i).lt.xtmin)goto 10 xplot = x(i) yplot = y(i) if(yplot.gt.ytmax)yplot = ytmax if(yplot.lt.ytmin)yplot = ytmin call movxy(xplot, yplot) call linxy(xplot, ytmin) 10 continue end c plot data as dashed lines subroutine odlplot(npts, x, y) include 'typecom' include 'robcom' double precision x(*), y(*) integer npts common/limits/xminp,xmaxp,yminp,ymaxp common/test/xtmin, xtmax, ytmin, ytmax ratio = 1.0d0 xscl = (0.821d0*fnplot/fmplot)**2 gs=(xmaxp-xminp)*(xmaxp-xminp)*xscl+(ymaxp-yminp)**2 do 10 i = 1,npts-1 x1 = x(i) x2 = x(i+1) y1 = y(i) y2 = y(i+1) if(y1.gt.ytmax)y1 = ytmax if(y1.lt.ytmin)y1 = ytmin if(x1.gt.xtmax)x1 = xtmax if(x1.lt.xtmin)x1 = xtmin if(y2.gt.ytmax)y2 = ytmax if(y2.lt.ytmin)y2 = ytmin if(x2.gt.xtmax)x2 = xtmax if(x2.lt.xtmin)x2 = xtmin fndash=sqrt(((x1-x2)*(x1-x2)*xscl+(y1-y2)*(y1-y2))/gs)*50.0d0 call dline(x1,y1,x2,y2,fndash,ratio) 10 continue end subroutine dlplot(npts, x, y) include 'typecom' include 'robcom' double precision x(*), y(*) common/limits/xminp,xmaxp,yminp,ymaxp common/test/xtmin, xtmax, ytmin, ytmax call dodash() do 8 i=1,npts if(x(i).ge.xtmin.and.x(i).le.xtmax.and.y(i).ge.ytmin + .and.y(i).le.ytmax)then call movxy(x(i),y(i)) goto 9 end if 8 continue 9 j = i do 10 i=j,npts xp = x(i) yp = y(i) if(xp.ge.xtmin.and.xp.le.xtmax.and.yp.ge.ytmin. +and.yp.le.ytmax)then call linxy(xp,yp) end if 10 continue call nodash() end c plot the data as lines joining points subroutine lplot(npts, x, y) include 'typecom' include 'robcom' double precision x(*), y(*) integer npts common/limits/xminp,xmaxp,yminp,ymaxp common/test/xtmin, xtmax, ytmin, ytmax c check on data limits do 8 i = 1, npts if(x(i) .ge. xtmin .and. x(i) .le. xtmax .and. y(i) .ge. ytmin + .and. y(i) .le. ytmax)then call movxy(x(i), y(i)) goto 9 end if 8 continue 9 j = i do 10 i = j, npts xp = x(i) yp = y(i) if(xp . ge. xtmin .and. xp. le. xtmax .and. yp .ge .ytmin + .and. yp .le. ytmax)then call linxy(xp, yp) end if 10 continue end c crude 3D plotting routines c only some of the standard 2D plot modes are supported c long term ambition is to have hidden line removal subroutine tdplot(npts, x, y, zfalse, +pmode, isymb, symb, sxl, syl, sxh, syh, ntickx, nticky, +paxes, ccode, labx, xlab, ylab, zlab, +tdfill) include 'typecom' include 'robcom' double precision x(*), y(*), zfalse(*) logical paxes, labx, ccode, tdfill common/limits/xminp,xmaxp,yminp,ymaxp common/limitz/zminp, zmaxp common/scales/xscale, yscale, zscale common/test/xtmin, xtmax, ytmin, ytmax character*(*) pmode, symb, xlab, ylab, zlab common/trig/cost1, sint1, cost2, sint2, cost3, sint3 c some local variables double precision x3min, x3max, y3min, y3max character*(60) buffer c degree to radian conversion c colours if needed c fudge for npz ntickz = (ntickx + nticky)/2 buffer = symb cost1 = cos(angle(1)*dtr) sint1 = sin(angle(1)*dtr) cost2 = cos(angle(2)*dtr) sint2 = sin(angle(2)*dtr) cost3 = cos(angle(3)*dtr) sint3 = sin(angle(3)*dtr) c find x/y/z minimum zmin = zfalse(1) zmax = zfalse(1) xhmin = x(1) yhmin = y(1) xhmax = x(1) yhmax = y(1) do 5 i = 2, npts zmin = min(zmin, zfalse(i)) zmax = max(zmax, zfalse(i)) xhmin = min(xhmin, x(i)) xhmax = max(xhmax, x(i)) yhmax = max(yhmax, y(i)) yhmin = min(yhmin, y(i)) 5 continue call linsc2(xhmin,xhmax,ntickx,xminp,xmaxp,npx,distx) call linsc2(yhmin,yhmax,nticky,yminp,ymaxp,npy,disty) call linsc2(zmin,zmax,ntickz,zminp,zmaxp,npz,distz) xscale = xmaxp - xminp yscale = ymaxp - yminp zscale = zmaxp - zminp c find max and min of x3 and y3 (3d equivalents of x and y) c and the corners of the box call rotate2(xminp, yminp, zminp, x3, y3) x3max = x3 x3min = x3 y3max = y3 y3min = y3 call rotate2(xmaxp, yminp, zminp, x3, y3) call maxmin3(x3max, x3min, y3max, y3min, x3, y3) call rotate2(xminp, ymaxp, zminp, x3, y3) call maxmin3(x3max, x3min, y3max, y3min, x3, y3) call rotate2(xminp, yminp, zmaxp, x3, y3) call maxmin3(x3max, x3min, y3max, y3min, x3, y3) call rotate2(xmaxp, ymaxp, zmaxp, x3, y3) call maxmin3(x3max, x3min, y3max, y3min, x3, y3) c the actual data do 10 i = 1, npts call rotate2(x(i), y(i), zfalse(i), x3, y3) call maxmin3(x3max, x3min, y3max, y3min, x3, y3) 10 continue c define plot limits sxh2=(sxh-sxl)*fiplot/fmplot+sxl sxl2=(sxh-sxl)*(fiplot-1.0d0)/fmplot+sxl syh2=(syh-syl)*fjplot/fnplot+syl syl2=(syh-syl)*(fjplot-1.0d0)/fnplot+syl call gap(sxl2, syl2, sxh2, syh2) if(tdfill)then call limsq(sxl2,syl2,sxh2,syh2, x3min, y3min, x3max, y3max) else call limsq(sxl2,syl2,sxh2,syh2, + -1.41d0, -1.41d0, 1.41d0, 1.41d0) end if c draw axes if desired if(paxes)then call rotate2(xminp, yminp, zminp, x3c, y3c) call movxy(x3c,y3c) call rotate2(xminp, yminp, zmaxp, x3, y3) call linxy(x3,y3) if(labx)then call rotate2(xminp, yminp, + (zminp+zmaxp)/2.0d0, x3, y3) call txtm(x3, y3, zlab, 60) end if call movxy(x3c,y3c) call rotate2(xminp, ymaxp, zminp, x3, y3) call linxy(x3,y3) if(labx)then call rotate2(xminp, (ymaxp+yminp)/2.0d0, + zminp, x3, y3) call txtm(x3, y3, ylab, 60) end if call movxy(x3c,y3c) call rotate2(xmaxp, yminp, zminp, x3, y3) call linxy(x3, y3) if(labx)then call rotate2((xminp+xmaxp)/2.0d0, + yminp, zminp, x3, y3) call txtm(x3, y3, xlab, 60) end if c complete the box call box3d(xminp, yminp, zminp, xmaxp, ymaxp, zmaxp) end if call rotate2(x(1), y(1), zfalse(1), x3, y3) call movxy(x3, y3) if(abs(rgbd(1)-rgb(1)).gt.2.or.abs(rgbd(2)-rgb(2)).gt.2.or. + abs(rgbd(3) -rgb(3)).gt.2)then call farkc(rgbd(1),rgbd(2),rgbd(3)) end if if(pmode.eq.'NICE')then do 50 i = 1, npts call rotate2(x(i), y(i), zfalse(i), x3, y3) if(CCODE)call datcol(zfalse(i), zminp, zscale) call symbol(x3, y3, isymb) 50 continue else if(pmode.eq.'CONTOUR')then call cont(x, y, zfalse, npts, ccode, .TRUE.) else if(pmode.eq.'SYMBOL')then do 55 i = 1, npts if(CCODE)call datcol(zfalse(i), zminp, zscale) call rotate2(x(i), y(i), zfalse(i), x3, y3) call ctext(x3,y3,buffer) 55 continue else if(pmode.eq.'PILLAR')then do 60 i = 1, npts if(CCODE)call datcol(zfalse(i), zminp, zscale) call rotate2(x(i), y(i), zfalse(i), x3, y3) call movxy(x3, y3) call rotate2(x(i), y(i), zminp, x3, y3) call linxy(x3, y3) 60 continue else if(pmode.eq.'LINES')then do 65 i = 1, npts if(CCODE)call datcol(zfalse(i), zminp, zscale) call rotate2(x(i), y(i), zfalse(i), x3, y3) call linxy(x3, y3) 65 continue else if(pmode.eq.'DASHEDLINES')then call dodash() do 70 i = 1, npts call rotate2(x(i), y(i), zfalse(i), x3, y3) call linxy(x3, y3) 70 continue call nodash() end if c make sure colour is back to normal call farkc(rgb(1),rgb(2),rgb(3)) end c 3D box subroutine box3d(x1, y1, z1, x2, y2, z2) include 'typecom' call mov3d(x1, y1, z1) call lin3d(x2, y1, z1) call lin3d(x2, y2, z1) call lin3d(x1, y2, z1) call lin3d(x1, y1, z1) call mov3d(x1, y1, z2) call lin3d(x2, y1, z2) call lin3d(x2, y2, z2) call lin3d(x1, y2, z2) call lin3d(x1, y1, z2) call mov3d(x1, y1, z1) call lin3d(x1, y1, z2) call mov3d(x1, y2, z1) call lin3d(x1, y2, z2) call mov3d(x2, y2, z1) call lin3d(x2, y2, z2) call mov3d(x2, y1, z1) call lin3d(x2, y1, z2) end subroutine mov3d(x, y, z) include 'typecom' call rotate2(x, y, z, x3, y3) call movxy(x3, y3) end subroutine lin3d(x, y, z) include 'typecom' call rotate2(x, y, z, x3, y3) call linxy(x3, y3) end subroutine rotate(x, y, z, cost1, sint1, cost2, sint2, + cost3, sint3, x3, y3) include 'typecom' c rotate about three axes in succession c y3 = z*cost1 - y*sint1 c x3 = x c z3 = y*cost1 + z*sint1 c next rotation c x3 = x3*cost2 + z3*sint2 c y3 = y3 c and last one c oldx3 = x3 c x3 = y3*sint3 + x3*cost3 c y3 = y3*cost3 - oldx3*sint3 double precision xt, yt, zt call rcex(xt, yt, zt, x, y, z) c try "standard" system (as defined in Foley and van Dam) call r1(xt, yt, zt, x3, y3, z3, cost1, sint1) call rcex(xt, yt, zt, x3, y3, z3) call r2(xt, yt, zt, x3, y3, z3, cost2, sint2) call rcex(xt, yt, zt, x3, y3, z3) call r3(xt, yt, zt, x3, y3, z3, cost3, sint3) end subroutine r1(xin, yin, zin, xout, yout, zout, c, s) include 'typecom' xout = xin*c - yin*s yout = xin*s + yin*c zout = zin end subroutine r2(xin, yin, zin, xout, yout, zout, c, s) include 'typecom' xout = xin yout = yin*c - zin*s zout = yin*s + zin*c end subroutine r3(xin, yin, zin, xout, yout, zout, c, s) include 'typecom' xout = xin*c - zin*s yout = yin zout = -xin*s + zin*c end subroutine rcex(x1, y1, z1, x2, y2, z2) include 'typecom' x1 = x2 y1 = y2 z1 = z2 end c scale to max and min subroutine rotate2(x, y, z, x3, y3) include 'typecom' common/trig/cost1, sint1, cost2, sint2, cost3, sint3 common/scales/xscale, yscale, zscale common/limits/xminp,xmaxp,yminp,ymaxp common/limitz/zminp, zmaxp x2 = (x-xminp)/xscale y2 = (y-yminp)/yscale z2 = (z-zminp)/zscale call rotate(x2, y2, z2, cost1, sint1, cost2, sint2, + cost3, sint3, x3, y3) end subroutine rotate3(x, y, z, x3, y3, threed) include 'typecom' logical threed include 'robcom' common/scales/xscale, yscale, zscale common/limits/xminp,xmaxp,yminp,ymaxp common/limitz/zminp, zmaxp if(threed)then cost1 = cos(angle(1)*dtr) sint1 = sin(angle(1)*dtr) cost2 = cos(angle(2)*dtr) sint2 = sin(angle(2)*dtr) cost3 = cos(angle(3)*dtr) sint3 = sin(angle(3)*dtr) x2 = (x-xminp)/xscale y2 = (y-yminp)/yscale z2 = (z-zminp)/zscale call rotate(x2, y2, z2, cost1, sint1, cost2, sint2, + cost3, sint3, x3, y3) else x3 = x y3 = y end if end c assign maximum and minimum values subroutine maxmin3(x3max, x3min, y3max, y3min, x3, y3) include 'typecom' x3max = max(x3max, x3) x3min = min(x3min, x3) y3max = max(y3max, y3) y3min = min(y3min, y3) end c interpolate colour value for data point c used to help show z scale in plot subroutine datcol(zdat, zminp, zscale) include 'typecom' include 'robcom' double precision zdat, zminp, zscale integer icol(3) double precision temp integer i do 10 i = 1, 3 temp = zdat-zminp temp = temp/zscale temp = temp*(rgbhi(i) - rgblo(i)) + rgblo(i) icol(i) = nint(temp) 10 continue call farkc(icol(1), icol(2), icol(3)) end c plot data as a histogram subroutine hplot(npts, x, y) include 'typecom' integer npts double precision x(*), y(*) include 'robcom' common/test/xtmin, xtmax, ytmin, ytmax common/limits/xminp,xmaxp,yminp,ymaxp double precision xp, yp C FIND START POINT INSIDE DATA LIMITS do 18 i=1,npts xp = x(i) yp = y(i) if(xp .ge. xtmin .and. xp .le. xtmax .and. yp .ge. ytmin + .and. yp .le. ytmax)then call movxy(x(i),y(i)) goto 19 end if 18 continue 19 J=I C EXTEND LINE BACKWARD IF POSSIBLE if(j.gt.1)xp=(x(i) + x(i-1))/2.0d0 if(j.eq.1)xp=x(i) - ((x(i+1)-x(i))/2.0d0) yp = y(i) if(xp .lt. xtmin)xp=xtmin if(yp .lt. ytmin)yp=ytmin if(xp .gt. xtmax)xp=xtmax if(yp .gt. ytmax)yp=ytmax call linxy(xp, yp) xp = x(j) yp = y(j) if(xp .lt. xtmin)xp = xtmin if(yp .lt. ytmin)yp = ytmin if(xp .gt. xtmax)xp = xtmax if(yp .gt. ytmax)yp = ytmax call movxy(xp, yp) do 110 i=j+1,npts xp=(x(i)+x(i-1))/2.0d0 yp=y(i-1) if(yp .gt. ytmax)yp=ytmax if(xp .lt. xtmin)xp=xtmin if(yp .lt. ytmin)yp=ytmin if(xp .gt. xtmax)goto 111 call linxy(xp,yp) yp=y(i) if(yp .gt. ytmax)yp=ytmax if(yp .lt. ytmin)yp=ytmin call linxy(xp,yp) xp=x(i) if(xp .gt. xtmax)goto 111 c for back to front data (e.g. after inverting x) if(xp .lt. xtmin) xp = xtmin call linxy(xp,yp) 110 continue C TERMINATE HISTOGRAM xp=x(npts)+((x(npts)-x(npts-1))/2.0d0) if(xp .gt. xtmax)xp=xtmax if(xp .lt. xtmin)xp = xtmin call linxy(xp,yp) goto 112 111 continue call linxy(xmaxp, yp) 112 continue end c plot data as a histogram c slightly different version that just uses error bars, c doesn't shrink or extend as "regular" version subroutine hplot2(npts, x, y, delx, dely) include 'typecom' integer npts double precision x(*), y(*), delx(*), dely(*) include 'robcom' common/test/xtmin, xtmax, ytmin, ytmax common/limits/xminp,xmaxp,yminp,ymaxp double precision xp, yp C FIND START POINT INSIDE DATA LIMITS do 18 i=1,npts xp = x(i) yp = y(i) if(xp .ge. xtmin .and. xp .le. xtmax .and. yp .ge. ytmin + .and. yp .le. ytmax)then call movxy(x(i),y(i)) goto 19 end if 18 continue 19 J=I C EXTEND LINE BACKWARD IF POSSIBLE xp = x(i) - delx(i) yp = y(i) if(xp .lt. xtmin)xp=xtmin if(yp .lt. ytmin)yp=ytmin if(xp .gt. xtmax)xp=xtmax if(yp .gt. ytmax)yp=ytmax call linxy(xp, yp) xp = x(i) + delx(i) yp = y(i) if(xp .lt. xtmin)xp = xtmin if(yp .lt. ytmin)yp = ytmin if(xp .gt. xtmax)xp = xtmax if(yp .gt. ytmax)yp = ytmax call linxy(xp, yp) do 110 i = j+1, npts xp=x(i) - delx(i) yp=y(i-1) if(yp .gt. ytmax)yp=ytmax if(xp .lt. xtmin)xp=xtmin if(yp .lt. ytmin)yp=ytmin if(xp .gt. xtmax)goto 111 call movxy(xp,yp) yp=y(i) if(yp .gt. ytmax)yp=ytmax if(yp .lt. ytmin)yp=ytmin call linxy(xp,yp) xp=x(i) if(xp .gt. xtmax)goto 111 call linxy(xp,yp) xp = xp + delx(i) if(xp .lt. xtmin)xp = xtmin if(yp .lt. ytmin)yp = ytmin if(xp .gt. xtmax)xp = xtmax if(yp .gt. ytmax)yp = ytmax call linxy(xp, yp) 110 continue C TERMINATE HISTOGRAM xp=x(npts)+delx(npts) if(xp .gt. xtmax)xp=xtmax if(xp .lt. xtmin)xp = xtmin call linxy(xp,yp) goto 112 111 continue c call linxy(xmaxp, yp) 112 continue end subroutine barplt(npts, xp, yp, xerr, yerr) include 'typecom' double precision xp(*), yp(*), xerr(*), yerr(*) include 'robcom' c do a 'bar' plot common/test/xtmin, xtmax, ytmin, ytmax do 1000 i = 1, npts x = xp(i) y = yp(i) xl = x - xerr(i) xr = x + xerr(i) yb = ytmin if(xl.gt.xtmax)goto 1000 if(xr.lt.xtmin)goto 1000 xl = MAX(xl, xtmin) xr = MIN(xr, xtmax) yt = MIN(ytmax, y) call fbox(xl, yb, xr, yt) 1000 continue end c make sure line width menus set correctly then call c actual routine subroutine lwid(i) include 'typecom' integer i call widchk(i) call lwidth(i) end c make sure font size menu items correctly set up then c call the real routine subroutine txtcall(size, chsz) include 'typecom' double precision size, chsz integer isize isize = nint(size) call txtchk(isize) call txtset(size) chsz = nint(size) end c zero a real array subroutine zain(ain, length) include 'typecom' double precision ain(*) integer i, length do 10 i = 1, length ain(i) = 0.0d0 10 continue end c zero an integer array subroutine zint(iin, length) include 'typecom' integer iin(*) integer i, length do 10 i = 1, length iin(i) = 0 10 continue end c bundles together two routines - FFT and Koji c Fourier routine. "itype" determines which one is used c "slow" Fourier transform c rips off an old routine Koji wrote c Returns "best period (bp) like c folding routines subroutine power(x, y, ndata, work1, work2, work3, yerr, + xout, yout, delx, dely, nout, + start, stop, nf, bp, + sxh, sxl, syh, syl, itype, + overwrite, plotft, ifft, iwtft) include 'typecom' double precision x(*), y(*), work1(*), work2(*), work3(*) double precision yerr(*) double precision xout(*), yout(*), delx(*), dely(*) c whether it's an FFT or FT integer itype c overwrite determines whether spectrum replaces data logical overwrite c plotft says whether we plot the spectrum or not logical plotft c What we do if data isn't 2**n points (0 = truncate, 1 = pad, c 2 = nearest of previous options c for adding "photons" double precision addup c for calculating max power/mean power double precision powersum c double precision templ integer ifft include 'robcom' include 'msizcom' c Only do transforms if there are at least 3 points if(ndata .le. 2)then call xtext('ERROR: Not enough data for transform') write(tstring, *)'Number of points = ', ndata call xtext(tstring) return end if fmax = 0.0d0 powersum = 0.0d0 fbest = fstart if(itype .eq. 1)then if(nf. gt. msiz)then call xtext('ERROR: TOO MANY FREQUENCIES') write(tstring, *)'Number requested = ', nf call xtext(tstring) write(tstring, *)' (Max allowed is ',msiz,')' call xtext(tstring) call xtext('Use ARRAYSIZE command to increase this') return end if if(start .ge. stop)then call xtext('ERROR: START VALUE > END VALUE') write(tstring, *)'( ', start, ',', stop, ' )' call xtext(tstring) return end if write(ifil8, *)'Periodogram:' write(ifil8, *)'start =',start,'stop =',stop write(ifil8, *)'Number of frequencies =',nf c convert periods to frequencies fstart = 1.0d0/stop fstop = 1.0d0/start c print*,'fstart = ',fstart, 'fstop = ', fstop fstep = (fstop - fstart)/nf c print*,'fstep = ', fstep write(tstring, *)'Calculating ', nf,' frequencies...' call xtext(tstring) if(iwtft .eq. 1)then call xtext('(Weighting data by errors)') end if fstart = fstart * twopi fstep = fstep * twopi call slwft(x, y, work1, work2, ndata, nf, fstart, + fstep, yerr, iwtft, work3) fstart = fstart/twopi fstep = fstep/twopi c crude normalization do 10 i = 1, nf work1(i) = work1(i)*work1(i) + work2(i)*work2(i) powersum = powersum + work1(i) if(work1(i) .gt. fmax)then fmax = work1(i) fbest = fstart + dble(i-1)*fstep c and period bp = 1.0d0/fbest end if c print*,work1(i) 10 continue c print*,'fmax =',fmax c FFT Part else if (itype .eq. 2)then c initial setting is complete data array n2use = ndata c put data into work array do 50 i = 1, n2use work1(i) = y(i) 50 continue c number of points must be 2**n ftest = n2use ftest = flog2(ftest) c this is floating point, do real comparison as integer itest = 2**nint(ftest) if((n2use - itest) .ne. 0)then call xtext('WARNING: Number of points not power of 2') write(tstring, *)'No. of points = ', n2use call xtext(tstring) c Do we truncate? Do this if 1) ifft = 0, or nearest, or arrays c aren't large enought to pad c Truncation mode (default) if(ifft .eq. 0)then n2use = 2**int(ftest) write(tstring, *)'truncating to ', n2use call xtext(tstring) c Nearest of truncate or pad (in log sense) else if(ifft .eq. 2)then if(itest .lt. n2use)then c we're going to truncate anyway n2use = 2**int(ftest) write(tstring, *)'truncating to ', n2use call xtext(tstring) else call fftpad(work1, ndata, n2use, itest, msiz, ftest) end if c Pad mode else if(ifft .eq. 1)then call fftpad(work1, ndata, n2use, itest, msiz, ftest) else c Unknown mode write(tstring, *)'Unknown FFTMODE of ', ifft call xtext(tstring) end if end if c crude check whether the data are approximately evenly spaced c need at least 3 data points if(n2use .ge. 3)then test = x(3) - x(2) - x(2) + x(1) test = test / (x(2) - x(1)) if(test . gt. 0.01d0)then call xtext('WARNING: DATA APPEAR TO BE UNEVENLY SPACED!') write(tstring, *)'(Test compared points 1, 2 and 3)' call xtext(tstring) end if c also check start and end times test = (ndata - 1)* (x(2) - x(1)) + x(1) c print*,'test = ',test, 'x(',n2use,') =', x(n2use) c print*, 'x(2) - x(1) =', x(2) - x(1) test = (test - x(ndata)) / (x(2) - x(1)) if (abs(test) .gt. 0.01d0)then call xtext('WARNING: DATA APPEAR TO BE UNEVENLY SPACED!') write(tstring,*)'(Test compared points 1, 2 and ', ndata,')' call xtext(tstring) end if end if nf = n2use/2 fstep = 1.0d0/(2.0d0*(nf)*(x(2) - x(1))) fstart = 1.0d0*fstep fstop = (nf )*fstep c print*,'fstep = ', fstep c number of "photons" in spectrum addup = 0.0d0 do 400 i = 1, n2use addup = addup + work1(i) 400 continue c print*,'addup = ',addup if(n2use .lt. 4)then call xtext('Not doing FFT because of small number of points') return end if call fft(work1, n2use, 1) work1(1) = 0.0d0 work1(2) = 0.0d0 c DONT USE THE LINES BELOW IF LEAHY NORMALIZATION IS USED c----- c was once creating an array of amplitudes, now change to power afactor = dble(nf)*dble(nf) afactor = 1.0d0/afactor do 60 i = 0, n2use, 2 temp = work1(i+1)*work1(i+1) temp = temp + work1(i+2)*work1(i+2) work1(i/2) = temp * afactor 60 continue c----- c Leahy normalization if more than one count c----- c Temporarily stop doing this c if (addup .gt. 1.0d0)then c c do 60 i = 0, n2use, 2 c temp = work1(i+1)*work1(i+1) c temp = temp + work1(i+2)*work1(i+2) c work1(i/2) = 2.0d0 * temp /addup c60 continue c c end if c----- c Leahy normalization do 61 i = 1, nf powersum = powersum + work1(i) if(work1(i) .gt. fmax)then fmax = work1(i) fbest = fstart + dble(i-1)*fstep c and period bp = 1.0d0/fbest c and bin number ibbest = i end if 61 continue else print*,'Unknown value of itype in Power = ', itype end if if(plotft)then sxh2=(sxh-sxl)*fiplot/fmplot+sxl sxl2=(sxh-sxl)*(fiplot-1.0d0)/fmplot+sxl syh2=(syh-syl)*fjplot/fnplot+syl syl2=(syh-syl)*(fjplot-1.0d0)/fnplot+syl c allow for window gaps call gap(sxl2, syl2, sxh2, syh2) call limit(sxl2,syl2,sxh2,syh2,fstart,0.0d0,fstop,fmax*1.2d0) call movxy(fstart, work1(1)) do 200 i = 1, nf value = fstart + (i-1)*fstep call linxy(value, work1(i)) 200 continue fpmax = 1.2d0 * fmax call boxm(fstart, 0.0d0, fstop, fpmax) write(tstring, *)'Peak at f = ', real(fbest), +'; Period = ',real(bp) call txtm(fstart, fpmax, tstring, nnl(tstring)) write(tstring, *)'Amp. = ',real(sqrt(fmax)) c Also get amplitude based on corrections for binning (frequency c dependent terms) and _average_ smearing of a signal c (caused by by a signal not corresponding precisely to a Fourier c frequency. call rtit(fstop, fpmax, tstring) write(tstring, *)real(fstart) call ctext(fstart, -fpmax/20.0d0, tstring) write(tstring, *)real(fstop) call ctext(fstop, -fpmax/20.0d0, tstring) write(tstring, *)'Frequency' call ctext(fstart + (fstop-fstart)/2.0d0, + -fpmax/10.0d0, tstring) c end of plot end if write(tstring, *) + 'Largest amplitude (amp^2) = ',sqrt(fmax), ' (',fmax,')' call xtext(tstring) call write8(tstring) c power ratio powersum = powersum/dble(nf) write(tstring, *) + 'Ratio of peak power to mean power =', fmax/powersum call xtext(tstring) call write8(tstring) c for FFT only if(itype .eq. 2)then templ = sqrt(2.0d0*(fmax -2.0d0)/addup) templ = templ*100.0d0 write(tstring, *)'"Leahy" power = ', templ, '%' call xtext(tstring) call write8(tstring) c add non-exact frequency effect templ = templ*sqrt(1.0d0/0.773d0) c add binning (frequency dependent) effect factor = dble(ibbest)/(dble(n2use)*pi) factor = sqrt(factor) * sin(factor) templ = templ*factor write(tstring, *)'With non-Fourier f correction =',templ call xtext(tstring) call write8(tstring) end if write(tstring, *)'Found at a period of ', bp call xtext(tstring) call write8(tstring) write(tstring, *)'equivalent to a frequency of ', fbest call xtext(tstring) call write8(tstring) c if we want to overwrite the data if(overwrite)then nout = nf do 100 i = 1, nf xout(i) = fstart + (i-1)*fstep yout(i) = work1(i) delx(i) = fstep/2.0d0 dely(i) = 0.0d0 100 continue end if end c From a routine by Koji Mukai subroutine slwft( X, Y, C, S, nData, nF, Fmin, dF, + yerr, iwtft, work) include 'typecom' C C Chebyshev recursion formulae used in Paul Murdin's version C accumulate numerical errors, and therefore now replaced C with a better routine taken from the Monro book. C C KM 11/7/85 C C Yerror array added - experiment with weighting of FT results c RC Dec. 1996 c Idea is to: c (1) Calculate weighted mean c (2) y -> (y - mean)/error**2 c (3) to avoid zero power problems subtract off newly calculated c mean of data set (bogus?) Integer nData, nF double precision X(*), Y(*), C(*), S(*), Fmin, dF double precision yerr(*) integer iwtft double precision Z0, dZ, cosK, sinK, dcos, dsin Integer J, K, Interval double precision work(*) Do K = 1, nF S(K) = 0.0d0 C(K) = 0.0d0 End Do c separate normalization into two bits to avoid (?) confusion if(iwtft .eq. 1)then call normftw(work, y, yerr, ndata) else if(iwtft .eq. 0)then call normft(work, y, ndata) else print*,'Unknown weighting = ',iwtft print*,'in slwft' end if c------ c------ c Print * If( nF .le. 100 ) Then Interval = 500 Else If( nF .le. 200 ) Then Interval = 500 ! 200 Else If( nF .le. 500 ) Then Interval = 500 ! 100 Else Interval = 500 ! 200 c were 1000, 500, 200 and 500 End If Do J = 1, nData Z0 = Fmin * X( J ) cosK = COS( Z0 ) sinK = SIN( Z0 ) dZ = dF * X( J ) dcos = COS( dZ ) dsin = SIN( dZ ) c C( 1 ) = C( 1 ) + cosK * Y( J ) c S( 1 ) = S( 1 ) + sinK * Y( J ) c--------- c new bit C( 1 ) = C( 1 ) + cosK * work(j) S( 1 ) = S( 1 ) + sinK * work(j) c end new bit c--------- Do K = 2, nF Call SWING( cosK, sinK, dcos, dsin ) c C( K ) = C( K ) + cosK * Y( J ) c S( K ) = S( K ) + sinK * Y( J ) c--------- c new bit C( K ) = C( K ) + cosK * work(j) S( K ) = S( K ) + sinK * work(j) c end new bit c--------- End Do c If ( MOD( J, Interval ) .eq. 0 ) Then c Write( *, 110 ) J, nData c110 Format( '+ Processed ', I5, ' out of ', I6, ' records' ) c End If End Do c Write( *, 120 ) nData c120 Format( '+ *Now finished processing', I7, ' records.' ) Z0 = 2.0d0 / dble( nData ) Do K = 1, nF C( K ) = C( K ) * Z0 S( K ) = S( K ) * Z0 End Do C End C C C subroutine SWING( cosA, sinA, cosB, sinB ) include 'typecom' C double precision half, one, one5 Parameter( half = 0.5d0, one = 1.0d0, one5 = 1.5d0 ) C double precision cosA, sinA, cosB, sinB C C cosA (out) = cos(A+B), sinA (out) = sin(A+B) C double precision D, Z1, Z2, T C D = cosB - one Z1 = cosA + D * cosA - sinA * sinB Z2 = sinA + D * sinA + cosA * sinB T = one5 - half * ( Z1 * Z1 + Z2 * Z2 ) cosA = T * Z1 sinA = T * Z2 C End c Normalize FT - weight by Errors subroutine normftw(work, y, yerr, ndata) include 'typecom' double precision work(*), y(*), yerr(*) integer ndata double precision mean, factor integer i double precision dmean, denoms, denom c prescale Y values by errors if required (RC) c For _all FTs_ we'll subtract off the mean dmean = 0.0d0 denoms = 0.0d0 do 10 i = 1, nData work(i) = y(i) if(yerr(i) .ne. 0.0d0)then denom = 1.0d0/(yerr(i)*yerr(i)) work(i) = y(i) * denom denoms = denoms + denom else work(i) = 0.0d0 end if dmean = dmean + work(i) 10 continue mean = dmean/denoms c For FT we subtract off mean of weighted values c instead of weighted mean (I think) c factor is to keep power spectrum appropriately weighted factor = dble(nData)/denoms do 11 i = 1, nData work(i) = (y(i) - mean)/(yerr(i)*yerr(i)) work(i) = factor * work(i) 11 continue c check what mean of work array is and subtract off c This stage can probably be removed dmean = 0.0d0 do 12 i = 1, ndata dmean = dmean + work(i) 12 continue mean = dmean/dble(ndata) do 13 i = 1, ndata work(i) = work(i) - mean 13 continue end c Normalize FT - no weighting by Errors subroutine normft(work, y, ndata) include 'typecom' double precision work(*), y(*) integer ndata double precision dmean double precision mean integer i dmean = 0.0d0 do 10 i = 1, nData work(i) = y(i) dmean = dmean + work(i) 10 continue mean = dmean/dble(ndata) c subtract it off do 13 i = 1, ndata work(i) = work(i) - mean 13 continue end c FFT routine based on algorithm in Numerical Recipes subroutine fft(yin, npts, isign) include 'typecom' double precision wr, wi, wpr, wpi, wtemp, theta double precision yin(*) n = npts / 2 c print*,'in fft, n = ', n c print*,'isign = ', isign theta = 6.28318530717959d0/2.0d0/dble(n) c1 = 0.5d0 if (isign .eq. 1) then c2 = -0.5d0 call fourr(yin, n, +1) else c2 = 0.5d0 theta = -theta endif wpr = -2.0d0*dsin(0.5d0*theta)**2 wpi = dsin(theta) wr = 1.0d0 + wpr wi = wpi n2p3 = 2 *n + 3 c do 11 i = 2, n/2+1 do 11 i = 2, n/2 i1 = 2 * i - 1 i2 = i1 + 1 i3 = n2p3 - i2 i4 = i3 + 1 wrs = sngl(wr) wis = sngl(wi) h1r = c1*(yin(i1) + yin(i3)) h1i = c1*(yin(i2) - yin(i4)) h2r = -c2*(yin(i2) + yin(i4)) h2i = c2*(yin(i1) - yin(i3)) yin(i1) = h1r + wrs*h2r - wis*h2i yin(i2) = h1i + wrs*h2i + wis*h2r yin(i3) = h1r - wrs*h2r + wis*h2i yin(i4) = -h1i + wrs*h2i + wis*h2r wtemp = wr wr = wr * wpr - wi * wpi +wr wi = wi * wpr + wtemp * wpi +wi 11 continue if (isign .eq. 1) then h1r = yin(1) yin(1) = h1r + yin(2) yin(2) = h1r - yin(2) else h1r = yin(1) yin(1) = c1*(h1r + yin(2)) yin(2) = c1*(h1r - yin(2)) call fourr(yin, n, -1) endif end subroutine fourr(data, nn, isign) include 'typecom' double precision wr, wi, wpr, wpi, wtemp, theta integer nn, isign double precision data(*) n = 2 * nn j = 1 do 11 i = 1, n, 2 if(j .gt. i)then tempr = data(j) tempi = data(j + 1) data(j) = data(i) data(j + 1) = data(i + 1) data(i) = tempr data(i + 1) = tempi endif m = n/2 1 if ((m .ge. 2) .and. (j .gt. m)) then j = j-m m = m/2 go to 1 endif j = j + m 11 continue mmax = 2 2 if (n .gt. mmax) then istep = 2 * mmax theta = 6.28318530717959d0/(isign*mmax) wpr = -2.d0 * dsin(0.5d0*theta)**2 wpi = dsin(theta) wr = 1.d0 wi = 0.d0 do 13 m = 1, mmax, 2 do 12 i = m, n, istep j = i + mmax c penny wise, pound foolish j1 = j + 1 i1 = i + 1 tempr = sngl(wr) * data(j) - sngl(wi) * data(j1) tempi = sngl(wr) * data(j1) + sngl(wi) * data(j) data(j) = data(i) - tempr data(j1) = data(i1) - tempi data(i) = data(i) + tempr data(i1) = data(i1) + tempi 12 continue wtemp = wr wr = wr * wpr - wi * wpi + wr wi = wi * wpr + wtemp * wpi + wi 13 continue mmax = istep go to 2 endif end c log base 2 of a number double precision function flog2(x) c inverse log base 10 of 2 double precision log102i double precision x data log102i/3.321928d0/ flog2 = log10(x) * log102i end c pad array with mean of data c assumes subroutine fftpad(work, ndata, n2use, iwant, msiz, ftest) include 'typecom' double precision work(*) integer ndata, n2use, iwant, msiz c ftest is log2 data size double precision ftest include 'robcom' double precision ave c see if our arrays are big enough if(iwant .gt. msiz)then call xtext('Arrays too small to pad') write(tstring, *)'Need a size of at least ',iwant call xtext(tstring) n2use = 2**int(ftest) write(tstring, *)'Truncating to ',n2use call xtext(tstring) return end if c get mean of data ave = 0.0d0 do 100 i = 1, ndata ave = ave + work(i) 100 continue ave = ave/dble(ndata) write(tstring, *)'Padding to ', iwant, ' points' call xtext(tstring) write(tstring, *)'With a value of ', ave call xtext(tstring) do 200 i = ndata+1, iwant work(i) = ave 200 continue n2use = iwant end