home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTH / 4THPROG.ZIP / PLTSPEC.FOR < prev    next >
Encoding:
Text File  |  1985-10-30  |  11.4 KB  |  334 lines
  1. $TITLE: 'PLOT AND MANIPULATE SPECTRUMS'
  2. $LARGE
  3. $NOFLOATCALLS
  4.       DIMENSION R(10,640),SCAT1(10,640),IPTS(10)
  5.       REAL DUMX(640),DUMY(640),KEVPC
  6.       CHARACTER*64 FNAME(10),ANSW(128)*1,CHR*1,HALDEV
  7.       INTEGER GMODE,GXMAX,GYMAX
  8.       INTEGER*4 IMAGE(8000),KEEP(8000)
  9.       LOGICAL IEXIST
  10. C
  11. C    This program was written to interactively process gamma-ray
  12. C    energy spectrums but is useful in any application where you 
  13. C    need to find the area under a peak, remove the background,
  14. C    fix lost/erroneous data points, or want to produce graphs
  15. C    with labels.
  16. C    The program is written in MS-FORTRAN and uses the MULTI-HALO
  17. C    graphics library extensively.
  18. C
  19. C  IPTS(j)         Number of points in spectrum j
  20. C  NSPT            Number of spectrums
  21. C  FNAME(j)        File name of spectrum j
  22. C  R,SCAT1         X and Y data for each spectrum
  23. C  DUMX,DUMY       Dummy vectors used when calling
  24. C                     some of the HALO routines.
  25. C  BASE,KEVPC      Two types of data files are used:
  26. C                  Type 1) file has only one column of 
  27. C                  numbers representing the Y values.  When file
  28. C                  type 1 is read the X values are computed assuming
  29. C                  the points are equally spaced.  BASE is the X value
  30. C                  of the first point and KEVPC is the distance between
  31. C                  successive points.
  32. C                  Type 2) file has two columns of numbers representing
  33. C                  the X and Y values.
  34. C                  The file type is indicated by placing an integer
  35. C                     1  or  2  on the first line of the file.
  36. C                  
  37. C
  38. C  IMAGE           Copy of the last graphics screen image
  39. C                      automatically saved each time.
  40. C  KEEP            Copy of a graphics screen image manually
  41. C                      saved by the user.
  42. C  GETCHR          Is a function which uses a machine language routine
  43. C                  to poll the keyboard for input.  The user does not
  44. C                  need to press the RETURN key as when a READ statement
  45. C                  is used.  Most important when using the graphics 
  46. C                  interactively----no extraneous characters are printed
  47. C                  on the screen when getting input from the keyboard.
  48. C
  49. C   LINKing--
  50. C       (.OBJ): PLTSPEC+EDITS+PLOT1+PLOT2+DOSFN+DOSFUNC+MANSPT+HALODVXX
  51. C(PLTSPEC.EXE): <RETURN>
  52. C    (NUL.MAP): <RETURN>
  53. C       (.LIB): A:+A:HALOF   (assuming libraries are on A:)
  54. C
  55. C  NOTE: Run-time disk needs the appropriate HALOxxx.DEV file.
  56. C
  57.       WRITE(0,'(A\)') ' Device driver name (ex. \d:HALOxxx.DEV\) = '
  58.       READ(0,'(A)') HALDEV
  59.       CALL SETDEV(HALDEV)
  60.       WRITE(0,9090)
  61.  9090 FORMAT(10X,'HALOIBM.DEV',/,25X,' 0:  320 x 200',/,
  62.      1       25X,' 1:  640 x 200',//,10X,'HALOSTB.DEV',/,
  63.      2       25X,' 0:  320 x 200 x 4',/,25X,' 1:  640 x 200 x 2',/,
  64.      3       25X,' 2:  320 x 200 x 16',/,25X,' 3:  640 x 200 x 4',/,
  65.      4       25X,'10:  640 x 352 x mono',/,1X)
  66.       WRITE(0,'(A\)') '  GRAPHICS MODE NUMBER: '
  67.       READ(0,*) GMODE
  68.       NSPT=0
  69.       DO 10 I=1,10
  70.       DO 5 J=1,640
  71.         R(I,J)=0.
  72.         SCAT1(I,J)=0.
  73.   5     CONTINUE
  74.  10     IPTS(I)=0
  75. C
  76.  9000 FORMAT(10X,'Read in a spectrum file',/,
  77.      1       10x,'Plot spectrums',/,
  78.      2       10x,'Set max/min limits',/,
  79.      4       10X,'Overlay Keep spectrum with another',/,
  80.      5       10X,'Last screen recall',/,
  81.      6       10X,'Graphics dump to printer',/,
  82.      7       10X,'Keep last screen',/,
  83.      8       10X,'Edit spectrum',/,
  84.      9       10X,'Manipulate last spectrum plotted',/,
  85.      3       10x,'Quit',/,1X)
  86.  20   WRITE(0,9000)
  87.       WRITE(0,'(A\)') '       COMMAND: '
  88.       CALL GETCHR(CHR)
  89.       ANSW(1)=CHR
  90.       WRITE(0,*)CHR
  91. C
  92. C
  93.       IF((ANSW(1).EQ.'R').OR.(ANSW(1).EQ.'r')) THEN
  94.    30   WRITE(0,'(A\)') ' ID NUMBER OF SPECTRUM TO READ '
  95.         READ(0,*,ERR=30) ID
  96.         WRITE(0,'(A\)') ' INPUT FILE NAME? '
  97.         READ(0,'(A)',ERR=30) FNAME(ID)
  98.         INQUIRE(FILE=FNAME(ID),EXIST=IEXIST)
  99.         IF(IEXIST) THEN
  100.           OPEN(7,FILE=FNAME(ID),STATUS='OLD')
  101.         ELSE
  102.           WRITE(0,*) ' FILE DOES NOT EXIST--- TRY AGAIN'
  103.           GO TO 30
  104.         ENDIF
  105. C--------------
  106.         READ(7,*) IFLTYP
  107.         IF(IFLTYP.EQ.1) THEN
  108. C--------------
  109.  31       WRITE(0,'(A\)')'  Energy of first channel and keV/channel: '
  110.           READ(0,*,ERR=31) BASE,KEVPC
  111. C
  112.           DO 90 J=1,640
  113.             READ(7,*,ERR=107,END=100) SCAT1(ID,J)
  114.             R(ID,J)=BASE + (J-1)*KEVPC
  115.   90      CONTINUE
  116.         ELSE
  117.           DO 95 J=1,640
  118.             READ(7,*,ERR=107,END=100) R(ID,J),SCAT1(ID,J)
  119.   95      CONTINUE
  120.         ENDIF
  121.  100    CLOSE(7,STATUS='KEEP')
  122.         J=J-1
  123.         IF(IPTS(ID).EQ.0) NSPT=NSPT + 1
  124.         IPTS(ID)=J
  125.         GO TO 20
  126. C
  127. C   error in reading file-- treat as if file was never opened.
  128. C
  129. 107     CLOSE(7,STATUS='KEEP')
  130.         IF(IPTS(ID).NE.0) NSPT=NSPT-1
  131.         IPTS(ID)=0
  132.         GO TO 20
  133.       ELSEIF((ANSW(1).EQ.'P').OR.(ANSW(1).EQ.'p')) THEN
  134.  249    DO 250 I=1,10
  135.           IF(IPTS(I).NE.0) WRITE(0,'(10X,I2,1X,A)') I,FNAME(I)
  136.  250    CONTINUE
  137.         WRITE(0,'(1X,/,1X,A\)') 'Spectrum number to plot '
  138.         READ(0,*,ERR=249) ID
  139. C
  140. C                          initialize graphics - 640x200 mono
  141.         CALL INITGR(GMODE)
  142.         CALL SETIEE(1)
  143.         CALL INQDRA(GXMAX,GYMAX)
  144. C
  145. C                           draw X- and Y-axis
  146.         CALL MOVABS(0,0)
  147.         CALL LNABS(0,GYMAX)
  148.         CALL LNABS(GXMAX,GYMAX)
  149.         JMAX=IPTS(ID)-1
  150.         CALL SETWOR(XMIN,RMIN,XMAX,RMAX)
  151. C
  152. C                           move graphics cursor to first point
  153. C
  154.         CALL MOVABS(R(ID,1),SCAT1(ID,1))
  155.         DO 255 J=1,IPTS(ID)
  156.           DUMY(J)=SCAT1(ID,J)
  157.  255      DUMX(J)=R(ID,J)
  158. C
  159. C                           graph the points
  160. C
  161.         CALL POLYLA(DUMX(2),DUMY(2),JMAX)
  162.         CALL SETGPR(2)
  163.         CALL WORLDO
  164. C
  165. C                           save a copy of the image
  166. C
  167.         CALL INQDRA(GXMAX,GYMAX)
  168.         CALL MOVEFR(0,0,GXMAX,GYMAX,IMAGE(1))
  169.         IDLSP=ID
  170.         CALL GETCHR(CHR)
  171.         CALL CLOSEG
  172.         GO TO 20
  173.       ELSEIF((ANSW(1).EQ.'S').OR.(ANSW(1).EQ.'s')) THEN
  174. C
  175. C     set upper and lower limits on X- and Y-axis
  176. C
  177.  1949   DO 1950 I=1,10
  178.           IF(IPTS(I).NE.0) WRITE(0,'(10X,I2,1X,A)') I,FNAME(I)
  179.  1950   CONTINUE
  180.         WRITE(0,'(10X,I2,1X,A)') 11,'To manually set limits'
  181.         WRITE(0,'(1X,/,1X,A\)') 'Spectrum number to get limits from '
  182.         READ(0,*,ERR=1949) ID
  183.         IF(ID.NE.11) THEN
  184.           JMAX=IPTS(ID)
  185.           XMAX=R(ID,1)
  186.           XMIN=R(ID,1)
  187.           RMAX=SCAT1(ID,1)
  188.           RMIN=SCAT1(ID,1)
  189.           DO 2000 I=1,JMAX
  190.             IF(SCAT1(ID,I).GT.RMAX) RMAX=SCAT1(ID,I)
  191.             IF(SCAT1(ID,I).LT.RMIN) RMIN=SCAT1(ID,I)
  192.             IF(R(ID,I).GT.XMAX) XMAX=R(ID,I)
  193.             IF(R(ID,I).LT.XMIN) XMIN=R(ID,I)
  194.  2000       CONTINUE
  195.           WRITE(0,*) ' Y-AXIS LIMITS ARE NOW ',RMAX,RMIN
  196.           WRITE(0,*) ' X-AXIS LIMITS ARE NOW', XMAX,XMIN
  197.         ELSE
  198.           WRITE(0,*) ' OLD Y-AXIS LIMITS ARE  ',RMAX,RMIN
  199.           WRITE(0,'(A\)') ' ENTER NEW LIMITS: '
  200.           READ(0,*,ERR=20) RMAX,RMIN
  201.           WRITE(0,*) ' OLD X-AXIS LIMITS ARE  ',XMAX,XMIN
  202.           WRITE(0,'(A\)') ' ENTER NEW LIMITS: '
  203.           READ(0,*,ERR=20) XMAX,XMIN
  204.         ENDIF
  205.         GO TO 20
  206.       ELSEIF((ANSW(1).EQ.'G').OR.(ANSW(1).EQ.'g')) THEN
  207. C
  208. C       graphics dump to an IDS Microprism printer
  209. C
  210.         WRITE(0,'(A\)') '  Half or Full height dump (H/F): '
  211.         READ(0,'(A)',ERR=20) ANSW(2)
  212.         WRITE(0,*) 'READY PRINTER AND PRESS RETURN'
  213.         READ(0,*)
  214.         CALL INITGR(GMODE)
  215.         CALL SETIEE(1)
  216.         CALL MOVETO(0,0,IMAGE(1),1)
  217.         IF((ANSW(2).EQ.'H').OR.(ANSW(2).EQ.'h')) THEN
  218. C
  219. C                    640x200 
  220.           CALL PLOT1
  221.         ELSE
  222. C
  223. C                    640x400 every other line is blank
  224.           CALL PLOT2
  225.         ENDIF
  226.         CALL GETCHR(CHR)
  227.         CALL CLOSEG
  228.         GO TO 20
  229.       ELSEIF((ANSW(1).EQ.'Q').OR.(ANSW(1).EQ.'q')) THEN
  230. C
  231. C                       quit- stop the program
  232.         STOP
  233.       ELSEIF((ANSW(1).EQ.'L').OR.(ANSW(1).EQ.'l')) THEN
  234. C
  235. C                           re-examine last screen image
  236.         CALL INITGR(GMODE)
  237.         CALL SETIEE(1)
  238.         CALL MOVETO(0,0,IMAGE(1),1)
  239. C
  240. C                          display image until a key is pressed
  241.         CALL GETCHR(CHR)
  242.         CALL CLOSEG
  243.         GO TO 20
  244.       ELSEIF((ANSW(1).EQ.'E').OR.(ANSW(1).EQ.'e')) THEN
  245. C
  246. C           edit ( change values of points) last spectrum plotted
  247. C
  248.         WRITE(0,*)'  EDITING SPECTRUM',IDLSP,'  ',FNAME(IDLSP)
  249.         CALL GETCHR(CHR)
  250.         CALL INITGR(GMODE)
  251.         CALL SETIEE(1)
  252.         CALL MOVETO(0,0,IMAGE(1),1)
  253.         CALL SETWOR(XMIN,RMIN,XMAX,RMAX)
  254. C
  255. C     interactive visual editing of the data points
  256. C
  257.         CALL EDITS(IDLSP,R,SCAT1,IPTS,IMAGE,GMODE)
  258. C
  259.         CALL WORLDO
  260.         CALL INQDRA(GXMAX,GYMAX)
  261.         CALL MOVEFR(0,0,GXMAX,GYMAX,IMAGE(1))
  262.         CALL CLOSEG
  263.         GO TO 20
  264.       ELSEIF((ANSW(1).EQ.'K').OR.(ANSW(1).EQ.'k')) THEN
  265. C
  266. C                           saves a copy of the last screen image
  267. C                           used with overlaying two spectrums.
  268.         CALL INITGR(GMODE)
  269.         CALL SETIEE(1)
  270.         CALL MOVETO(0,0,IMAGE(1),1)
  271.         CALL INQDRA(GXMAX,GYMAX)
  272.         CALL MOVEFR(0,0,GXMAX,GYMAX,KEEP(1))
  273.         CALL CLOSEG
  274.         GO TO 20
  275.       ELSEIF((ANSW(1).EQ.'O').OR.(ANSW(1).EQ.'o')) THEN
  276. C
  277. C                           overlay one spectrum on another
  278.  549    DO 550 I=1,10
  279.           IF(IPTS(I).NE.0) WRITE(0,'(10X,I2,1X,A)') I,FNAME(I)
  280.  550    CONTINUE
  281.         WRITE(0,'(1X,/,1X,A\)') 'Spectrum number to plot '
  282.         READ(0,*,ERR=549) ID
  283.         CALL INITGR(GMODE)
  284.         CALL SETIEE(1)
  285.         CALL MOVETO(0,0,KEEP(1),1)
  286. C
  287. C  NOTE: the X and Y axis limits are set to those associated with
  288. C        the last 'Set limits' command which may not be the same
  289. C        as those used to draw the graph saved in the "KEEP" array
  290. C        or those of the spectrum being overlayed.
  291. C
  292.         CALL SETWOR(XMIN,RMIN,XMAX,RMAX)
  293.         CALL MOVABS(R(ID,1),SCAT1(ID,1))
  294.         JMAX=IPTS(ID)-1
  295.         DO 555 J=1,IPTS(ID)
  296.           DUMY(J)=SCAT1(ID,J)
  297.  555      DUMX(J)=R(ID,J)
  298.         CALL POLYLA(DUMX(2),DUMY(2),JMAX)
  299.         CALL SETGPR(2)
  300.         CALL WORLDO
  301.         CALL INQDRA(GXMAX,GYMAX)
  302.         CALL MOVEFR(0,0,GXMAX,GYMAX,IMAGE(1))
  303.         CALL GETCHR(CHR)
  304.         CALL CLOSEG
  305.         GO TO 20
  306.       ELSEIF((ANSW(1).EQ.'M').OR.(ANSW(1).EQ.'m')) THEN
  307.         WRITE(0,*)'  MANIPULATING SPECTRUM',IDLSP,'  ',FNAME(IDLSP)
  308.         DO 600 I=1,1000
  309.  600     RJER=8.8*9./5.
  310.         CALL INITGR(GMODE)
  311.         CALL SETIEE(1)
  312.         CALL MOVETO(0,0,IMAGE(1),1)
  313. C
  314.         CALL SETWOR(XMIN,RMIN,XMAX,RMAX)
  315. C
  316. C       MANSPT allows placing labels on a graph, fitting a second 
  317. C              order curve to the spectrum background, removing the
  318. C              background, calculating the counts under a peak, saving
  319. C              a copy of the screen image on the printer, and 
  320. C              saving the background corrected data in a file.
  321. C
  322.         CALL MANSPT(IDLSP,R,SCAT1,IPTS,IMAGE,GMODE)
  323. C
  324.         CALL WORLDO
  325.         CALL INQDRA(GXMAX,GYMAX)
  326.         CALL MOVEFR(0,0,GXMAX,GYMAX,IMAGE(1))
  327.         CALL CLOSEG
  328.         GO TO 20
  329.       ELSE
  330.         CALL CLOSEG
  331.         GO TO 20
  332.       ENDIF
  333.       END
  334.