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