sw/rescue/lab8e_goettingen/disk2_11/rkb/paroff/haescm.ft

   1 C       -+-+-+-+-+ \ e  H A E S C M . F T  \ e -+-+-+-+-+
   2 C
   3 C       COMPUTE SOME VALUES IN THE X Y AREA
   4 C
   5 C       UP  2-MAR-84  ADDS ANGULAR DISTRBUTION HISTOGRAM
   6 C       UP  4-MAR-84  PLOTS THE ANGULAR DISTRIBUTION HISTOGRAM
   7 C
   8         SUBROUTINE SCMPUT
   9         INCLUDE HAEBUF.FI
  10         INCLUDE HAETTY.FI
  11         INCLUDE HAEGSA.FI
  12         INCLUDE HAEPTI.FI
  13         INCLUDE HAPPEN.FI
  14         INCLUDE HAECSZ.FI
  15         INCLUDE HABRK.FI
  16         REAL XYSAM,TIME,ATAN2,COND
  17         INTEGER TOASCI,NCOND
  18         EXTERNAL XYSAM,TIME,TOASCI,NCOND,ATAN2,COND
  19 C
  20 \f       REAL DELX,DELY,XY
  21         REAL XAMED,YAMED,               @ ARITHMETISCHER MITTELWERT ALLER PUNKTE
  22      *  CPUTIM,                         @ HERE WE STORE THE COMPUTERS TIME NEEDED
  23      *  MERADI,DMERAD,                  @ MITTLERE GEOMETRISCHE ENTFERNUNG (XAMED,YAMED) ZU ALLEN PUNKTEN
  24 C                                       @ DAS IST EIN KREIS UM DEN MITTELPUNKT MIT RADIUS MERADI,
  25 C                                       @ DMERAD IST DIE STANDARTABWEICHUNG DAZU
  26      *  MEXY,                           @ MITTLERER WEG ZWISCHEN ZWEI MESSPUNKTEN
  27      *  DEMEXY,                         @ STANDARTABWEICHUNG DAZU
  28      *  PERIOD,                         @ 1/SAMRAT
  29      *  FLACH                           @ MERADI**2 * 3.1416 :  DIE FLAECHE DES KREISES UM XAMED,YAMED
  30      *  ,ANGHIS(36),WINK(36)            @ ANGULAR DISTRIBUTION HISTOGRAM
  31         INTEGER I,J,POS1,POSIT,
  32      *  S0,S1,                          @ BEGIN, END OF EXAMINATION
  33      *  MM,SS,MMM,SSS,TE,               @ SOME TEMPS
  34      *  CA,CO                           @ CHANNEL OF ABZISSE, ORDINATE
  35 C
  36 C       STATEMENT FUNCTIONS:
  37 C
  38         XSAM(IP1)=XYSAM(IP1,CA)
  39         YSAM(IP1)=XYSAM(IP1,CO)
  40         ANGLE(IP3)=MOD(INT((ABS(ATAN2(YSAM(IP3)-YAMED,XSAM(IP3)-XAMED))
  41      *   +COND(YSAM(IP3)-YAMED.LT.0,3.1416))*5.7297+.5),36)+1 @ THIS REALLY WORKS
  42         PHYTAG(P1,P2)=SQRT(P1*P1+P2*P2)
  43         DIST (IP3)   =PHYTAG(XAMED-XSAM(IP3),YAMED-YSAM(IP3))
  44         XCON(IP4)=COND(ABS(XSAM(IP4)-XSAM(IP4-1)).GT.3,
  45      *   XSAM(IP4)-XSAM(IP4-1))
  46         YCON(IP4)=COND(ABS(YSAM(IP4)-YSAM(IP4-1)).GT.3,
  47      *   YSAM(IP4)-YSAM(IP4-1))
  48 C
  49         IF (SAMCNT.EQ.0) RETURN         @ NO SAMPLE TOKEN, RETURN
  50         IF (OPTION.NE.12) CALL ASKHIM(4) @ IF SC* THEN DO NOT ASK FOR BEGIN,SPAN
  51         IF (BREAK(11)) GOTO 70          @ USER GETS RID OF PROGRAM
  52         CPUTIM=TIME(I)                  @ GET THE TIME OF START
  53 C
  54         WRITE (3,2)
  55         CA=0
  56         CO=1
  57         DO 50 POS1=BEGIN*SAMRAT,ENDS*SAMRAT-1,OVRLAP*SAMRAT
  58         POSIT=POS1+SPAN*SAMRAT          @ THIS IS THE END OF THE ACTUAL SCAN
  59         IF (POSIT.GT.ENDS*SAMRAT) GOTO 60 @ SKIP LAST OVERLAPPING COMPUTATION
  60         S0=POS1
  61         S1=S0+SPAN*SAMRAT
  62         XAMED=XSAM(S0)
  63         YAMED=YSAM(S0)
  64 C
  65         DO 10 I=S0+1,S1,1
  66         XAMED=XSAM(I)+XAMED
  67         YAMED=YSAM(I)+YAMED
  68 10      CONTINUE
  69         YAMED=YAMED/(S1-S0)             @ ARITHMETISCHER MITTELWERT
  70         XAMED=XAMED/(S1-S0)
  71 C
  72 C       BERECHNE DIE WEGLAENGE
  73 C
  74         XY=0
  75         IF(BREAK(11)) GOTO 70
  76         DO 15 I=S0+1,S1                 @ S1-S0 PUNKTE
  77
  78 15      XY=XY+PHYTAG(XCON(I),YCON(I))
  79 C
  80 C       MITTLERE WEGLAENGE IM EINEM MESSINTERVALL (1/SAMRAT SEC.)
  81 C
  82         MEXY=XY/(S1-S0)                 @ MITTLERE WEGLAENGE ZWISCHEN ZWEI PUNKTEN --> MEXY
  83         DEMEXY=0                        @ UND JETZT BERECHNEN WIR DIE STANDARTABWEICHUNG
  84         DO 17 I=S0+1,S1                 @ S1-S0 PUNKTE
  85 17      DEMEXY=(MEXY-PHYTAG(XCON(I),YCON(I)))**2+DEMEXY @ S1-S0-1 PUNKTE NUR!
  86         DEMEXY=SQRT(DEMEXY/(S1-S0-1))   @ STANDARTABWEICHUNG --> DEMEXY
  87 C
  88 C       MITTLERER KREIS UM (XAMED,YAMED)
  89 C
  90         MERADI=0
  91         DO 20 I=S0,S1                   @ S1-S0+1 PUNKTE
  92 20      MERADI=MERADI+DIST(I)
  93         MERADI=MERADI/(1+S1-S0)         @ MITTLERER KREIS
  94 C
  95 C       STANDARTABWEICHUNG ZU MERADI
  96 C
  97         DMERAD=0
  98         DO 30 I=S0,S1                   @ S1-S0+1 PKT
  99 30      DMERAD=(MERADI-DIST(I))**2+DMERAD
 100         DMERAD=SQRT(DMERAD/(S1-S0))
 101 C
 102 C       ANGULAR DISTRIBUTION HISTOGRAM ( HUFSCHMIDT ET AL.)
 103         IF (BREAK(11)) GOTO 70
 104         CALL MOVE (-36,0,ANGHIS)
 105         DO 35 I=S0,S1                   @ ANOTHER SCAN
 106         TE=ANGLE(I)
 107 35      ANGHIS(TE)=ANGHIS(TE)+1         @ HERE WE COMPUTE THE HISTOGRAM
 108 \fC
 109 C
 110 C       PREPARE THE OUTPUT LINE
 111 C
 112         MM=TOASCI(S0/SAMRAT/60)
 113         SS=TOASCI(MOD(S0/SAMRAT,60))
 114         MMM=TOASCI((S1/SAMRAT)/60)
 115         SSS=TOASCI(MOD(S1/SAMRAT,60))
 116 C
 117         PERIOD=1./SAMRAT
 118         FLACH =MERADI**2*3.141593
 119         WRITE (3,1) LABEL,MM,SS,MMM,SSS,
 120      *  XAMED,YAMED,MERADI,FLACH,DMERAD,XY,PERIOD,MEXY,DEMEXY
 121         WRITE (3,4) ((ANGHIS(I+J),J=1,28,9),I=0,8)
 122 C
 123 C       HERE WE PLOT THE SWAY DIRECTION HISTOGRAM
 124         DO 110 I=1,36
 125 110     WINK(I)=(I*10.-5.)*.017453      @ FILL 0 TO 2*PI IN 5 DEGREE STEPS
 126         IF (BREAK(11)) GOTO 70
 127         CALL TWAIT (10000)              @ WAIT FOR THE PRINTER TO COMPLETE
 128         CALL STPLT
 129         CALL XYPLOT (XOFSET,YOFSET,-PENUP)
 130         CALL XYPLOT (XLEN/2.,(YLEN-1.)/2.,-PENUP) @ SET ORIGIN FOR POLAR PLOT RTN
 131         CALL POLAR (ANGHIS,WINK,36,1,0,3,0,(YLEN-1.)/200.)
 132         CALL XYPLOT(-XLEN/2.,-(YLEN-1.)/2.,-PENUP)
 133         CALL LABPLT(0,-YOFSET,ENDS)
 134         IF (BREAK(11)) GOTO 120
 135         CALL SYMBOL (-YZ,2*YZ,YZ,'SWAY DIRECTION HISTOGRAM',90.,25)
 136         CALL SYMBOL(YZ,YLEN+YZ,YZ,LABEL,0,MIN0(42,LABCNT*6)) @ PLOT 42 CHARACTERS INTO THE FIRST LINE OF LABEL
 137         IF (LABCNT.GT.7)
 138      *  CALL SYMBOL (YZ,YLEN,YZ,LABEL(8),0,LABCNT*6-42) @ AND THE REMAINING INTO THE NEXT LINE
 139 120     CALL XYPLOT (38.,25.,PENUP)
 140         CALL EXPLT
 141         IF (BREAK(11)) GOTO 70          @ USER GET'S RID OF PLOTTER
 142 C
 143 C
 144 50      CONTINUE
 145 60      CONTINUE
 146         CPUTIM=TIME(I)-CPUTIM           @ COMPUTE EXECUTION TIME
 147         WRITE (TTO,3) CPUTIM
 148 70      RETURN
 149 1       FORMAT (1X,10A6,T100,A2,1H:,A2,' --> ',A2,1H:,A2,
 150      *  /'  MITTE X',F7.1,'  Y',F7.1,5X,
 151      *  'RADIUS',F6.1,' FLAECHE',F8.0,' S.A.',F5.1,
 152      *  3X,'WEG',F7.0,3X,'FUER',F6.3,' SEC.  WEG',F6.3,' S.A.',F6.2)
 153 2       FORMAT (/50(2H -)/1X,13A6)
 154 3       FORMAT (' RECHENZEIT:',F6.1,' SEC')
 155 4       FORMAT ('0 HISTOGRAMM DER WINKELVERTEILUNG:  ARCTAN (YI/XI)',
 156      *   '  POSITION HISTOGRAM',//,22X,
 157      *   ,' 1. QUADRANT   2. QUADRANT   3. QUADRANT   4. QUADRANT '//
 158      *  ,2(5(21X,4F12.0/)/))
 159         END
 160 \1a