home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTH / 4THPROG.ZIP / MANSPT.FOR < prev    next >
Encoding:
Text File  |  1985-10-30  |  14.1 KB  |  450 lines
  1. $TITLE:' MANIPUTLATE AND LABEL SPECTRUM'
  2. $LARGE R,SPECT1,IPTS,IMAGE
  3. $NOFLOATCALLS
  4.       SUBROUTINE MANSPT(IDLSP,R,SPECT1,IPTS,IMAGE,GMODE)
  5.       REAL*4 R(10,640),SPECT1(10,640)
  6.       CHARACTER*1 RCHR,ANSW,F2CHR*30,LABEL*60,FONTXT*20,CHR
  7.       LOGICAL IEXIST
  8.       INTEGER IPTS(10),GMODE,GXM,GYM
  9.       INTEGER*4 IMAGE(8000),SAVE(8000)
  10.       INTEGER*2  CH,L1,R1,FH,ST,BFIT,EFIT
  11.        REAL DUMSP(640)
  12.        REAL*8 CDBLP(3,3),RSDBLP(3),QDBLP,QQDBLP
  13. C
  14. C    The following functions can be performed by this subroutine:
  15. C
  16. C  *fit the spectrum background to a 2nd order curve
  17. C
  18. C   commands-   B   begin fit at current cursor position
  19. C               L   set current cursor position as the left
  20. C                     side of the peak.
  21. C               R   set current cursor position as the right
  22. C                     side of the peak.
  23. C               E   end fit at current cursor position.
  24. C               C   compute the polynomial fit.
  25. C               X   remove the background from the data between
  26. C                     the points B and E and display the
  27. C                     corrected data.
  28. C               O   overlay the background curve between the 
  29. C                     points B and E.
  30. C               A   compute the area between the points L and R
  31. C                     taking into account the background curve
  32. C                     fit.  (Function X does not have to have been
  33. C                     done for this to work, however, function C
  34. C                     must have been performed.)
  35. C
  36. C     The data points between B and E (except those points between
  37. C     L and R) are fitted to a 2nd order polynomial
  38. C               y= AA*X**2 + BB*X + CC
  39. C     using a weighted least squares routine.  The weighting is 
  40. C     just 1/(square root of the counts).  NOTE: the fit will fail
  41. C     unless B is less than L is less than R is less than E.
  42. C
  43. C *start over with the original data
  44. C       command-    '-'  pressing this key erases everything and
  45. C                           starts over.
  46. C
  47. C *exit the subroutine
  48. C     commands-    Q   quits and returns to main program as if 
  49. C                        nothing had been done to the data.
  50. C                  F   finish manipulating the data and display-
  51. C                      returns to the main program with the screen
  52. C                      image intact and with any changes made to the
  53. C                      data (e.g. background corrected, etc.).
  54. C
  55. C *obtain printer output
  56. C     command-     G   dump screen image to printer
  57. C
  58. C
  59. C *write the data to a disk file
  60. C         command-     W   you will be prompted for a file name.
  61. C
  62. C *place labels on the screen image
  63. C    commands-
  64. C             T   first key to push.  you will be prompted
  65. C                 for the text of the label.  NOTE the first and
  66. C                 last charachters must be unique:
  67. C
  68. C                     \peak area = 19654 counts\
  69. C                     !the exclamation points are unique!
  70. C
  71. C
  72. C             8   move up        | With the NumLck key
  73. C             9   move up fast   | is lit on the IBM
  74. C             6   move right     | keyboard this key
  75. C             +   move right fast| arrangement really
  76. C             2   move down      | does make some
  77. C             3   move down fast | sense
  78. C             4   move left      |
  79. C        <RETURN> move left fast |
  80. C
  81. C             L   display the label starting at cursor position.
  82. C             E   erase the label. If the label is not displayed
  83. C                   the results are unpredictable.
  84. C             Z   allows you to change the height, width, writing
  85. C                   direction(path), and writing mode.
  86. C                   height and width can be 1,2,3 but not 1.5,2.3
  87. C                   direction   0   horizontal
  88. C                               1   at 90 deg to horizontal
  89. C                               2   upside down and backwards
  90. C                               3   at 270 deg to horizontal
  91. C                   mode 0 or 1     unboxed or boxed.
  92. C             Q   quit the labeling routine.
  93. C
  94.  10    L1=1
  95.        R1=10
  96.        CH=1
  97.        BFIT=1
  98.        EFIT=IPTS(IDLSP)
  99.        IINT=IPTS(IDLSP)
  100. C
  101. C               default label text parameters
  102.       IHT=1
  103.       IWDTH=1
  104.       IPTH=0
  105.       MODE=0
  106. C
  107.        DO 50 J=1,IINT
  108.  50      DUMSP(J)=SPECT1(IDLSP,J)
  109. C
  110.        CALL INQWOR(XW1,YW1,XW2,YW2)
  111.        CALL WORLDO
  112.        CALL INQDRA(GXM,GYM)
  113.        CALL MOVEFR(0,0,GXM,GYM,SAVE(1))
  114.        CALL INITHC(10,10,1)
  115.        CALL SETWOR(XW1,YW1,XW2,YW2)
  116.  100   CALL INKEY(RCHR)
  117.        X=R(IDLSP,CH)
  118.        Y=DUMSP(CH)
  119.        CALL MOVHCA(X,Y)
  120.        IF(RCHR.EQ.CHAR(0)) GO TO 100
  121.        IF(RCHR.EQ.'4') THEN
  122. C
  123. C        CURSOR LEFT
  124. C
  125.          CH=CH-1
  126.          IF (CH.LT.1) CH=1                                    
  127.          GO TO 100
  128.        ENDIF
  129.        IF(RCHR.EQ.CHAR(13)) THEN
  130.          CH=CH-5
  131.          IF(CH.LT.1) CH=1
  132.          GO TO 100
  133.        ENDIF
  134.        IF(RCHR.EQ.'6') THEN
  135. C
  136. C        CURSOR RIGHT
  137. C
  138.          CH=CH + 1
  139.          IF (CH.GT.IINT) CH=IINT
  140.          GO TO 100
  141.        ENDIF
  142.        IF(RCHR.EQ.'+') THEN
  143.          CH=CH+5
  144.          IF(CH.GT.IINT) CH=IINT
  145.          GO TO 100
  146.        ENDIF
  147.        IF((RCHR.EQ.'R').OR.(RCHR.EQ.'r')) R1=CH
  148.        IF((RCHR.EQ.'L').OR.(RCHR.EQ.'l')) L1=CH
  149.        IF((RCHR.EQ.'B').OR.(RCHR.EQ.'b')) BFIT=CH
  150.        IF((RCHR.EQ.'E').OR.(RCHR.EQ.'e')) EFIT=CH
  151.        IF((RCHR.EQ.'C').OR.(RCHR.EQ.'c')) GO TO 450
  152.        IF((RCHR.EQ.'X').OR.(RCHR.EQ.'x')) GO TO 1130
  153.        IF(RCHR.EQ.'-') GO TO 200
  154.        IF((RCHR.EQ.'A').OR.(RCHR.EQ.'a')) GO TO 1210
  155.        IF((RCHR.EQ.'O').OR.(RCHR.EQ.'o')) GO TO 1300
  156.        IF((RCHR.EQ.'W').OR.(RCHR.EQ.'w')) GO TO 1380
  157.        IF((RCHR.EQ.'T').OR.(RCHR.EQ.'t')) GO TO 1600
  158.        IF((RCHR.EQ.'G').OR.(RCHR.EQ.'g')) GO TO 1700
  159.        IF((RCHR.EQ.'F').OR.(RCHR.EQ.'f')) GO TO 1900
  160.        IF((RCHR.EQ.'Q').OR.(RCHR.EQ.'q')) GO TO 2000
  161.        GO TO 100
  162. C
  163. C    START OVER
  164. C
  165.  200   CALL WORLDO
  166.        CALL MOVETO(0,0,SAVE(1),1)
  167.        CALL INITHC(10,10,1)
  168.        GO TO 10
  169. C
  170. C       DO 2ND ORDER FIT
  171.  450   CONTINUE
  172.        DO 490 J=1,3
  173.          RSDBLP(J)=0.
  174.          DO 490 I=1,3
  175.          CDBLP(J,I)=0.
  176.  490   CONTINUE
  177.        FH=L1
  178.        ST=BFIT
  179.        IF (ST.GT.FH) ST=1
  180. C
  181. C       REM COMPUTE COEF. AND RHS
  182. C
  183.        DO 500 J=ST,FH
  184.          QDBLP=0.
  185.          IF(DUMSP(J).NE.0.)QDBLP=1./SQRT(DUMSP(J))
  186.          CDBLP(3,3)=CDBLP(3,3)+QDBLP
  187.          RSDBLP(3)=RSDBLP(3) +DUMSP(J)*QDBLP
  188.          QDBLP=QDBLP*R(IDLSP,J)
  189.          CDBLP(2,3)=CDBLP(2,3) + QDBLP
  190.          RSDBLP(2)=RSDBLP(2) + DUMSP(J)*QDBLP
  191.          QDBLP=QDBLP*R(IDLSP,J)
  192.          CDBLP(1,3)=CDBLP(1,3) + QDBLP
  193.          RSDBLP(1)=RSDBLP(1) + DUMSP(J)*QDBLP
  194.          QDBLP=R(IDLSP,J)*QDBLP
  195.          CDBLP(1,2)=CDBLP(1,2) + QDBLP
  196.          CDBLP(1,1)=CDBLP(1,1) + QDBLP*R(IDLSP,J)
  197.  500   CONTINUE
  198. C
  199.        FH=EFIT
  200.        ST=R1
  201.        IF (FH.LT.ST) FH=IINT
  202. C
  203.        DO 600 J=ST,FH
  204.          QDBLP=0.
  205.          IF(DUMSP(J).NE.0.) QDBLP=1./SQRT(DUMSP(J))
  206.          CDBLP(3,3)=CDBLP(3,3)+QDBLP
  207.          RSDBLP(3)=RSDBLP(3) +DUMSP(J)*QDBLP
  208.          QDBLP=QDBLP*R(IDLSP,J)
  209.          CDBLP(2,3)=CDBLP(2,3) + QDBLP
  210.          RSDBLP(2)=RSDBLP(2) + DUMSP(J)*QDBLP
  211.          QDBLP=QDBLP*R(IDLSP,J)
  212.          CDBLP(1,3)=CDBLP(1,3) + QDBLP
  213.          RSDBLP(1)=RSDBLP(1) + DUMSP(J)*QDBLP
  214.          QDBLP=R(IDLSP,J)*QDBLP
  215.          CDBLP(1,2)=CDBLP(1,2) + QDBLP
  216.          CDBLP(1,1)=CDBLP(1,1) + QDBLP*R(IDLSP,J)
  217.  600   CONTINUE
  218. C
  219.        CDBLP(2,1)=CDBLP(1,2)
  220.        CDBLP(2,2)=CDBLP(1,3)
  221.        CDBLP(3,1)=CDBLP(1,3)
  222.        CDBLP(3,2)=CDBLP(2,3)
  223.        DO 740 J=1,2
  224.          QQDBLP=CDBLP(J,J)                                          
  225.          DO 660 I=J,3                                           
  226.            CDBLP(J,I)=CDBLP(J,I)/QQDBLP                                  
  227.  660     CONTINUE
  228.          RSDBLP(J)=RSDBLP(J)/QQDBLP                                    
  229.          DO 740 I=J+1,3                                       
  230.            QQDBLP=CDBLP(I,J)                                          
  231.            DO 720 M=1,3
  232.              CDBLP(I,M)=CDBLP(I,M) - QQDBLP*CDBLP(J,M)                        
  233.  720       CONTINUE
  234.            RSDBLP(I)=RSDBLP(I) - QQDBLP*RSDBLP(J)                           
  235.  740   CONTINUE
  236.        CC=RSDBLP(3)/CDBLP(3,3)                                    
  237.        BB=RSDBLP(2) - CC*CDBLP(2,3)                               
  238.        AA=RSDBLP(1) - BB*CDBLP(1,2) - CC*CDBLP(1,3)
  239. C
  240.        CALL INQWOR(XW1,YW1,XW2,YW2)
  241.        CALL WORLDO
  242.        CALL INQDRA(GXM,GYM)
  243.        CALL MOVEFR(0,0,GXM,GYM,IMAGE(1))
  244.        WRITE(0,*) AA,BB,CC
  245.        CALL GETCHR(CHR)
  246.        CALL MOVETO(0,0,IMAGE(1),1)
  247.        CALL SETWOR(XW1,YW1,XW2,YW2)
  248.        GO TO 100
  249. C
  250. C       REM SUBTRACT BACKGROUND
  251. C
  252.  1130  CONTINUE
  253.        DO 1190 J=BFIT,EFIT
  254.        DUMSP(J)=DUMSP(J) - ((AA*R(IDLSP,J) + BB)*R(IDLSP,J) +CC)
  255.        IF (DUMSP(J).LT.0) DUMSP(J)=0
  256.  1190  CONTINUE
  257.        CALL REPLOT(IDLSP,R,DUMSP,IPTS,GMODE)
  258.         CALL WORLDO
  259.         CALL INQDRA(GXM,GYM)
  260.         CALL MOVEFR(0,0,GXM,GYM,IMAGE(1))
  261.        CALL SETWOR(XW1,YW1,XW2,YW2)
  262.        X=R(IDLSP,CH)
  263.        Y=DUMSP(CH)
  264.        CALL MOVHCA(X,Y)
  265.        GO TO 100
  266. C
  267. C       REM COMPUTE AREA OF PEAK
  268. C
  269.  1210  AR=0
  270.        DO 1250 J=L1,R1
  271.          AR=SPECT1(IDLSP,J) - ((AA*R(IDLSP,J) +BB)*R(IDLSP,J) + CC)+AR
  272.  1250  CONTINUE
  273. C
  274.        CALL INQWOR(XW1,YW1,XW2,YW2)
  275.        CALL WORLDO
  276.        CALL INQDRA(GXM,GYM)
  277.        CALL MOVEFR(0,0,GXM,GYM,IMAGE(1))
  278.        WRITE(0,*)'PEAK AREA =',AR
  279.        CALL GETCHR(CHR)
  280.        CALL MOVETO(0,0,IMAGE(1),1)
  281.        CALL SETWOR(XW1,YW1,XW2,YW2)
  282.        GO TO 100
  283. C
  284. C       REM OVERLAY BACKGROUND
  285. C
  286.  1300  CONTINUE
  287.        CALL SETLNS(2)
  288.        DO 1360 J=BFIT,EFIT
  289.          Y=(AA*R(IDLSP,J) +BB)*R(IDLSP,J) + CC
  290.          X=R(IDLSP,J)
  291.          IF(J.EQ.BFIT) THEN
  292.            CALL MOVABS(X,Y)
  293.          ELSE
  294.            CALL LNABS(X,Y)
  295.          ENDIF
  296.  1360  CONTINUE
  297.        CALL SETLNS(1)
  298.        GO TO 100
  299. C
  300. C       REM WRITE RESULTS
  301. C
  302.  1380  CONTINUE
  303.        CALL INQWOR(XW1,YW1,XW2,YW2)
  304.        CALL WORLDO
  305.        CALL INQDRA(GXM,GYM)
  306.        CALL MOVEFR(0,0,GXM,GYM,IMAGE(1))
  307.        WRITE(0,'(A\)')'  OUTPUT FILENAME: '
  308.        READ(0,'(A)') F2CHR
  309.        INQUIRE(FILE=F2CHR,EXIST=IEXIST)
  310.        IF(IEXIST) THEN
  311.          WRITE(0,'(A\)') ' File exists- overwrite it (Y/N): '
  312.          READ(0,'(A)')ANSW
  313.          IF((ANSW.EQ.'Y').OR.(ANSW.EQ.'y')) THEN
  314.            OPEN(4,FILE=F2CHR,STATUS='OLD')
  315.          ELSE
  316.            GO TO 1380
  317.          ENDIF
  318.        ELSE
  319.          OPEN(4,FILE=F2CHR,STATUS='NEW')
  320.        ENDIF
  321. C
  322.  1383  WRITE(0,'(A\)') ' 1 OR 2 COLUMN TYPE FILE (1/2): '
  323.        READ(0,*,ERR=1383) IFLTYP
  324.        IF(IFLTYP.EQ.1) THEN
  325.          WRITE(4,*,ERR=1435) IFLTYP
  326.          DO 1430 J=1,IINT
  327.            WRITE(4,'(1X,F8.1)',ERR=1435) DUMSP(J)
  328.  1430    CONTINUE
  329.        ELSEIF(IFLTYP.EQ.2) THEN
  330.          WRITE(4,*,ERR=1435) IFLTYP
  331.          DO 1432 J=1,IINT
  332.            WRITE(4,'(1X,F10.5,1X,F8.1)',ERR=1435) R(IDLSP,J),DUMSP(J)
  333.  1432    CONTINUE
  334.        ELSE
  335.          GO TO 1383
  336.        ENDIF
  337.  1435  CLOSE(4,STATUS='KEEP')
  338.        CALL MOVETO(0,0,IMAGE(1),1)
  339.        CALL SETWOR(XW1,YW1,XW2,YW2)
  340.        GO TO 100
  341. C
  342. C   LABELING ROUTINE
  343. C
  344.  1600 CONTINUE
  345.       CALL WORLDO
  346.       CALL INQDRA(GXM,GYM)
  347.       CALL MOVEFR(0,0,GXM,GYM,IMAGE(1))
  348.  1605 WRITE(0,'(A\)') ' TYPE IN TEXT FOR LABEL:'
  349.       READ(0,'(A)') LABEL
  350.       CALL MOVETO(0,0,IMAGE(1),1)
  351.       CALL SETTEX(IHT,IWDTH,IPTH,MODE)
  352.       CALL SETTCL(1,0)
  353.       CALL INITTC(8,8,1)
  354.       CALL SETXOR(1)
  355.  1610 CALL GETCHR(RCHR)
  356.       IF(RCHR.EQ.'2') THEN
  357.          CALL MOVTCR(0,1)
  358.       ELSEIF(RCHR.EQ.'3') THEN
  359.          CALL MOVTCR(0,5)
  360.       ELSEIF(RCHR.EQ.'6') THEN
  361.          CALL MOVTCR(1,0)
  362.       ELSEIF(RCHR.EQ.'+') THEN
  363.          CALL MOVTCR(5,0)
  364.       ELSEIF(RCHR.EQ.'8') THEN
  365.          CALL MOVTCR(0,-1)
  366.       ELSEIF(RCHR.EQ.'9') THEN
  367.          CALL MOVTCR(0,-5)
  368.       ELSEIF(RCHR.EQ.'4') THEN
  369.          CALL MOVTCR(-1,0)
  370.       ELSEIF(RCHR.EQ.CHAR(13)) THEN
  371.          CALL MOVTCR(-5,0)
  372.       ELSEIF((RCHR.EQ.'L').OR.(RCHR.EQ.'l')) THEN
  373.          CALL INQTCU(IX,IY,ICOL)
  374.          CALL TEXT(LABEL)
  375.       ELSEIF((RCHR.EQ.'E').OR.(RCHR.EQ.'e')) THEN
  376.          CALL MOVTCA(IX,IY)
  377.          CALL TEXT(LABEL)
  378.          CALL MOVTCA(IX,IY)
  379.       ELSEIF((RCHR.EQ.'Z').OR.(RCHR.EQ.'z')) THEN
  380.          CALL INQDRA(GXM,GYM)
  381.          CALL MOVEFR(0,0,GXM,GYM,IMAGE(1))
  382.  1606   WRITE(0,'(A\)') ' HEIGHT(int),WIDTH(int),PATH(int),MODE(int):'
  383.          READ(0,*,ERR=1606) IHT,IWDTH,IPTH,MODE
  384.          CALL MOVETO(0,0,IMAGE(1),1)
  385.          CALL SETTEX(IHT,IWDTH,IPTH,MODE)
  386.       ELSEIF((RCHR.EQ.'Q').OR.(RCHR.EQ.'q')) THEN
  387.          GO TO 1620
  388.       ENDIF
  389.       GO TO 1610
  390.  1620 CONTINUE
  391.       CALL DELTCU
  392.       CALL SETXOR(0)
  393.       CALL SETWOR(XW1,YW1,XW2,YW2)
  394.       GO TO 100
  395. C**********************  END LABELING ROUTINE
  396. C
  397. C   GRAPH TO PRINTER
  398. C
  399.  1700 CONTINUE
  400.       CALL WORLDO
  401.       CALL DELHCU
  402.        CALL INQDRA(GXM,GYM)
  403.        CALL MOVEFR(0,0,GXM,GYM,IMAGE(1))
  404.        WRITE(0,*) ' READY PRINTER AND PRESS ANY KEY'
  405.        CALL GETCHR(RCHR)
  406.        WRITE(0,'(A\)')' Half of Full height dump (H/F): '
  407.        READ(0,'(A)') RCHR
  408.        CALL MOVETO(0,0,IMAGE(1),1)
  409.       IF((RCHR.EQ.'H').OR.(RCHR.EQ.'h')) THEN
  410.         CALL PLOT1
  411.       ELSE
  412.         CALL PLOT2
  413.       ENDIF
  414.       CALL INITHC(10,10,1)
  415.       CALL SETWOR(XW1,YW1,XW2,YW2)
  416.       GO TO 100
  417. C
  418. C   EXIT
  419. C
  420.  1900  CONTINUE
  421.        DO 1905 J=1,IPTS(IDLSP)
  422.  1905    SPECT1(IDLSP,J)=DUMSP(J)
  423. C
  424. C   QUIT
  425. C
  426.  2000  CONTINUE
  427.        RETURN
  428.        END
  429.        SUBROUTINE REPLOT(ID,R,DUMSP,IPTS,GMODE)
  430.         INTEGER GMODE
  431.         DIMENSION R(10,640),DUMSP(640),IPTS(10),DUMX(640)
  432.         CALL INQWOR(XW1,YW1,XW2,YW2)
  433.         CALL WORLDO
  434.         CALL CLOSEG
  435.         CALL INITGR(GMODE)
  436.         CALL SETIEE(1)
  437.         CALL INQDRA(IXM,IYM)
  438.         CALL MOVABS(0,0)
  439.         CALL LNABS(0,IYM)
  440.         CALL LNABS(IXM,IYM)
  441.         JMAX=IPTS(ID)-1
  442.         CALL SETWOR(XW1,YW1,XW2,YW2)
  443.         CALL MOVABS(R(ID,1),DUMSP(1))
  444.         DO 255 J=1,IPTS(ID)
  445.  255      DUMX(J)=R(ID,J)
  446.         CALL POLYLA(DUMX(2),DUMSP(2),JMAX)
  447.         CALL SETGPR(2)
  448.         RETURN
  449.         END
  450.