BASIC HdB

10 REM Probleem van Labrujere
20 REM ======================
50 REM MESH GENERATOR
60 REM --------------
130 INPUT "Aantal rasterlijnen in radiale richting (7) : ",NDR
140 INPUT "Aantal rasterlijnen in omtreks-richting (5) : ",NPH
150 IF NDR<=0 THEN NDR=7
160 IF NPH<=0 THEN NPH=5
170 REM
180 REM ... Initialisatie ...
190 LIMC=NDR*NPH : REM aantal rasterpunten
193 LIMT=(NDR-1)*(NPH-1) : REM aantal elementen
197 LIMV=NDR*(NPH-1)+NPH*(NDR-1) : REM aantal knooppunten
200 DIM NOV(LIMV,2) : DIM MAV(LIMC,2) : REM Snelheden
210 LIMB=2*(NDR-1)+2*(NPH-1) : REM aantal randpunten
220 REM
230 PRINT "Algemene gegevens ..."
240 OPEN "R",#1,"HDBGENER.BIN",20
250 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$, 4 AS F5$
260 LSET F1$=MKI$(LIMC) : LSET F2$=MKI$(LIMT)
265 LSET F3$=MKI$(LIMV)
270 LSET F4$=MKI$(LIMB) : LSET F5$=MKI$(NPH)
280 PUT #1,1 : CLOSE #1
290 REM
300 PRINT "Coordinaten ..."
310 OPEN "R",#1,"HDBCOORD.BIN",8
320 FIELD #1, 4 AS X$, 4 AS Y$
330 PI=4*ATN(1) : DPHI=.5*PI/(NPH-1)
340 R=10 : P=R^(1/(NDR-1))-1 : REM Laatste radiale meshlijn
350 FOR I=1 TO NDR
360 R=(1+P)^(I-1) : REM Leidt tot equivalente elementen
370 FOR J=1 TO NPH
380 PHI=(J-1)*DPHI : K=(I-1)*NPH+J
390 X=R*COS(PHI) : Y=R*SIN(PHI)
400 LSET X$=MKS$(X) : LSET Y$=MKS$(Y)
410 PUT #1,K : NEXT J : NEXT I
420 CLOSE #1
430 REM
440 PRINT "Topologie ..."
450 OPEN "R",#1,"HDBTOPOL.BIN",16
460 FIELD #1, 4 AS F1$, 4 AS F2Y$, 4 AS F3$, 4 AS F4$
470 FOR I=1 TO NDR-1 : FOR J=1 TO NPH-1
480 K=(I-1)*(NPH-1)+J
490 KON2=(I-1)*NPH+J : KON1=(I-1)*NPH+J+1
500 KON3=I*NPH+J : KON4=I*NPH+J+1
510 REM PRINT KON1,KON2,KON3,KON4
520 LSET F1$=MKI$(KON1) : LSET F2$=MKI$(KON2)
530 LSET F3$=MKI$(KON3) : LSET F4$=MKI$(KON4)
540 PUT #1,K : NEXT J : NEXT I
550 CLOSE #1
560 REM
570 PRINT "Onbekenden ..."
580 FOR I=1 TO NDR : FOR J=1 TO NPH-1
590 NR1=(I-1)*NPH+J : NR2=(I-1)*NPH+J+1
600 K=(2*NPH-1)*(I-1)+J : NOV(K,1)=NR1 : NOV(K,2)=NR2
610 REM PRINT K,NOV(K,1),NOV(K,2)
620 NEXT J : NEXT I
630 FOR I=1 TO NDR-1 : FOR J=1 TO NPH
640 NR1=(I-1)*NPH+J : NR2=I*NPH+J
650 K=(2*NPH-1)*I+J-NPH : NOV(K,1)=NR1 : NOV(K,2)=NR2
660 REM PRINT K,NOV(K,1),NOV(K,2)
670 NEXT J : NEXT I
680 REM
690 FOR K=1 TO LIMV
700 I=NOV(K,1) : J=NOV(K,2)
710 II=I : IF J<I THEN II=J
720 JJ=ABS(I-J) : IF JJ>2 THEN JJ=2
730 MAV(II,JJ)=K : NEXT K
740 REM
750 OPEN "R",#1,"HDBUNKNO.BIN",8
760 FIELD #1, 4 AS F1$, 4 AS F2$
770 FOR I=1 TO LIMV
780 LSET F1$=MKI$(NOV(I,1)) : LSET F2$=MKI$(NOV(I,2))
790 PUT #1,I : NEXT I
800 FOR I=1 TO LIMC
810 LSET F1$=MKI$(MAV(I,1)) : LSET F2$=MKI$(MAV(I,2))
820 PUT #1,I+LIMV : NEXT I
830 CLOSE #1
840 PRINT "Randvoorwaarden ..."
850 OPEN "R",#1,"HDBBOUND.BIN",12
860 FIELD #1, 4 AS F1$, 4 AS F2Y$, 4 AS F3$
870 FOR I=1 TO NDR-1 : REM ... Instroming ...
880 NOB1=I*NPH : NOB2=(I+1)*NPH : NOB3=1
890 REM PRINT NOB1,NOB2,NOB3
900 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
905 LSET F3$=MKI$(NOB3)
910 PUT #1,I : NEXT I : L=NDR-1
920 FOR I=1 TO NPH-1 : REM ... Cirkel ...
930 NOB1=I : NOB2=I+1 : NOB3=2
940 REM PRINT NOB1,NOB2,NOB3
950 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
955 LSET F3$=MKI$(NOB3)
960 PUT #1,I+L : NEXT I : L=L+NPH-1
970 FOR I=1 TO NDR-1 : REM ... Symmetrie ...
980 NOB1=NPH*(I-1)+1 : NOB2=NPH*I+1 : NOB3=1
990 REM PRINT NOB1,NOB2,NOB3
1000 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
1005 LSET F3$=MKI$(NOB3)
1010 PUT #1,I+L : NEXT I : L=L+NDR-1
1020 FOR I=1 TO NPH-1 : REM ... Uitstroming ...
1030 NOB1=NPH*(NDR-1)+I : NOB2=NPH*(NDR-1)+I+1 : NOB3=3
1040 REM PRINT NOB1,NOB2,NOB3
1050 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
1055 LSET F3$=MKI$(NOB3)
1060 PUT #1,I+L : NEXT I : L=L+NPH-1 : CLOSE #1
1070 END

10 REM Probleem van Labrujere
20 REM ======================
30 REM HOOFDPROGRAMMA
40 REM --------------
45 PRINT "Gegevens inlezen ..."
100 OPEN "R",#1,"HDBGENER.BIN",20
110 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$, 4 AS F5$
120 GET #1,1 : LIMC=CVI(F1$) : LIMT=CVI(F2$)
125 LIMV=CVI(F3$) : LIMB=CVI(F4$)
130 NPH=CVI(F5$)
140 CLOSE #1 : NN=2*LIMV : NB1=4*NPH
150 REM
160 DIM X(LIMC),Y(LIMC),MAV(LIMC,2),S(NN,NB1),B(NN)
170 REM
180 REM Coordinaten:
190 OPEN "R",#1,"HDBCOORD.BIN",8
200 FIELD #1, 4 AS F1$, 4 AS F2$
210 FOR I=1 TO LIMC : GET #1,I
220 X(I)=CVS(F1$) : Y(I)=CVS(F2$)
230 NEXT I : CLOSE #1
240 REM
250 REM Nummering onbekenden:
260 OPEN "R",#1,"HDBUNKNO.BIN",8
270 FIELD #1, 4 AS F1$, 4 AS F2$
280 FOR I=1 TO LIMC : GET #1,I+LIMV
290 MAV(I,1)=CVI(F1$) : MAV(I,2)=CVI(F2$)
300 NEXT I : CLOSE #1
305 REM
320 REM Grote matrix en vector schooonmaken:
330 FOR I=1 TO NN : FOR J=1 TO NB1 : S(I,J)=0
335 NEXT J : B(I)=0 : NEXT I
340 REM
350 DIM NO(4),E(8,8),R(2),NR(8),A(8)
490 OPEN "R",#1,"HDBTOPOL.BIN",16
500 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$ : GET #1,1
510 NO(1)=CVI(F1$) : NO(2)=CVI(F2$)
515 NO(3)=CVI(F3$) : NO(4)=CVI(F4$)
520 PRINT "Repeteerbaar element ..." : GOSUB 1060
530 REM
540 PRINT "Assemblage procedure voor";LIMT;"elementen ..."
550 FOR F=1 TO LIMT : GET #1,F : PRINT F;
560 NO(1)=CVI(F1$) : NO(2)=CVI(F2$)
565 NO(3)=CVI(F3$) : NO(4)=CVI(F4$)
570 REM
580 REM Snelheids-nummering:
590 FOR K=1 TO 4 : L=(K MOD 4)+1
600 II=NO(K) : IF NO(L)<NO(K) THEN II=NO(L)
610 JJ=ABS(NO(K)-NO(L)) : IF JJ>2 THEN JJ=2
620 IJ=MAV(II,JJ) : NR(2*K-1)=2*IJ-1 : NR(2*K)=2*IJ
630 NEXT K
640 REM Element intellen:
650 FOR I=1 TO 8 : IG=NR(I)
660 FOR J=I TO 8 : JG=NR(J)
670 II=IG : IF JG<IG THEN II=JG
680 JJ=ABS(IG-JG)+1
690 S(II,JJ)=S(II,JJ)+E(I,J)
700 NEXT J : NEXT I
710 NEXT F : CLOSE #1 : PRINT : REM Einde bulk.
720 REM
730 PRINT "Randvoorwaarden ..."
740 OPEN "R",#1,"HDBBOUND.BIN",12
750 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$
760 FOR F=1 TO LIMB : GET #1,F
770 MU=CVI(F1$) : NU=CVI(F2$) : KIND=CVI(F3$)
780 GOSUB 1970
790 REM Snelheids-nummering:
800 II=MU : IF NU<MU THEN II=NU
810 JJ=ABS(MU-NU) : IF JJ>2 THEN JJ=2
820 IJ=MAV(II,JJ) : NR(1)=2*IJ-1 : NR(2)=2*IJ
830 FOR I=1 TO 2 : IG=NR(I)
840 FOR J=I TO 2 : JG=NR(J)
850 II=IG : IF JG<IG THEN II=JG
860 JJ=ABS(IG-JG)+1
870 S(II,JJ)=S(II,JJ)+E(I,J)
880 NEXT J
890 B(IG)=B(IG)+R(I)
900 NEXT I
910 NEXT F : CLOSE #1 : REM Einde totale assemblage.
920 REM
930 GOSUB 2180 : REM Los de vergelijkingen op ...
940 REM
950 OPEN "R",#1,"HDBVELOC.BIN",8
960 FIELD #1, 4 AS F1$, 4 AS F2$
970 FOR K=1 TO LIMV
980 LSET F1$=MKS$(B(2*K-1)) : LSET F2$=MKS$(B(2*K))
990 REM LPRINT B(2*K-1),B(2*K)
1000 PUT #1,K : NEXT K : CLOSE #1
1010 END
1020 REM =========================
1030 REM  Eindige Elementen Matrix
1040 REM =========================
1050 REM Maak schoon:
1060 FOR I=1 TO 8 : FOR J=1 TO 8 : E(I,J)=0 : NEXT J : NEXT I
1570 REM Lokaliseer gegevens:
1580 X1=X(NO(1)) : X2=X(NO(2)) : X3=X(NO(3)) : X4=X(NO(4))
1590 Y1=Y(NO(1)) : Y2=Y(NO(2)) : Y3=Y(NO(3)) : Y4=Y(NO(4))
1600 X21=X2-X1 : X32=X3-X2 : X43=X4-X3 : X14=X1-X4
1610 Y21=Y2-Y1 : Y32=Y3-Y2 : Y43=Y4-Y3 : Y14=Y1-Y4
1620 REM Onsamendrukbaar:
1630 A(1)=-Y21 : A(3)=-Y32 : A(5)=-Y43 : A(7)=-Y14
1640 A(2)=+X21 : A(4)=+X32 : A(6)=+X43 : A(8)=+X14
1650 GOSUB 1900 : REM Voeg toe aan element ...
1660 REM Rotatievrij:
1670 A(1)=+X21 : A(3)=+X32 : A(5)=+X43 : A(7)=+X14
1680 A(2)=+Y21 : A(4)=+Y32 : A(6)=+Y43 : A(8)=+Y14
1690 GOSUB 1900 : REM Voeg toe aan element ...
1700 REM Isoparametriek:
1710 A(1)=+1 : A(2)=0 : A(3)=-1 : A(4)=0
1715 A(5)=+1 : A(6)=0 : A(7)=-1 : A(8)=0
1720 GOSUB 1900 : REM Voeg toe aan element ...
1730 A(1)=0 : A(2)=+1 : A(3)=0 : A(4)=-1
1735 A(5)=0 : A(6)=+1 : A(7)=0 : A(8)=-1
1740 GOSUB 1900 : REM Voeg toe aan element ...
1750 RETURN
1760 REM --------------------
1770 REM Voeg toe aan element
1780 REM --------------------
1790 REM Spoorwegen:
1900 SUM=0 : FOR I=1 TO 8 : SUM=SUM+A(I)*A(I)
1905 NEXT I : SUM=1/SUM
1905 REM Kleinste kwadraten methode:
1910 FOR I=1 TO 8 : FOR J=1 TO 8
1920 E(I,J)=E(I,J)+A(I)*A(J)*SUM
1930 NEXT J : NEXT I
1940 RETURN
1950 REM ===============
1960 REM Randvoorwaarden
1970 REM ===============
1980 FOR I=1 TO 2 : R(I)=0 : FOR J=1 TO 2 : E(I,J)=0
1985 NEXT J : NEXT I
1990 REM
2000 ON KIND GOTO 2030,2060,2120
2010 REM
2020 REM ... V = 0 ...
2030 E(2,2)=1 : RETURN
2040 REM
2050 REM ... Y.U - X.V = 0 ...
2060 X21=X(NU)-X(MU) : Y21=Y(NU)-Y(MU)
2070 SUM=X21*X21+Y21*Y21
2080 E(1,1)=Y21*Y21/SUM : E(2,2)=X21*X21/SUM
2090 E(1,2)=-X21*Y21/SUM : E(2,1)=-X21*Y21/SUM
2100 RETURN
2110 REM ... U = 1 ...
2120 E(1,1)=1 : R(1)=1 : RETURN
2130 REM
2140 REM ==============================================
2150 REM  OPLOSSEN VAN LINEAIRE VERGELIJKINGEN VOOR
2160 REM  POSITIEF DEFINIETE SYMMETRISCHE BAND MATRICES
2170 REM ==============================================
2180 IF NN<2 THEN PRINT "GAUSS: NN=";NN : END
2190 IF NB1<2 THEN PRINT "GAUSS: NB1=";NB1 : END
2200 REM
2210 PRINT "Vegen voor";NN;"vergelijkingen ..."
2220 II=NN-1 : KK=NB1-1
2230 FOR I=1 TO II : DIAG=S(I,1) : PRINT I;
2240 IF DIAG=0 THEN PRINT "Nul op de diagonaal" : END
2250 IF I>NN-NB1 THEN KK=NN-I
2260 FOR K=1 TO KK
2270 PIVOT=S(I,K+1)/DIAG : IF PIVOT=0 THEN GOTO 2300
2280 B(I+K)=B(I+K)-PIVOT*B(I) : JJ=KK-K+1
2290 FOR J=1 TO JJ : S(I+K,J)=S(I+K,J)-PIVOT*S(I,K+J) : NEXT J
2300 NEXT K : NEXT I
2310 PRINT "Terugsubstitutie ..."
2320 KK=NB1-1 : FOR N=2 TO NN : I=NN-N+2 : PRINT I;
2330 B(I)=B(I)/S(I,1) : IF I<NB1 THEN KK=I-1
2340 FOR K=1 TO KK : B(I-K)=B(I-K)-B(I)*S(I-K,K+1) : NEXT K
2350 NEXT N : B(1)=B(1)/S(1,1) : RETURN
2360 REM

10 REM Probleem van Labrujere
20 REM ======================
50 REM Plot mesh & snelheden
60 REM ---------------------
100 REM Algemeen ...
110 OPEN "R",#1,"HDBGENER.BIN",20
120 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$, 4 AS F5$
130 GET #1,1 : LIMC=CVI(F1$) : LIMT=CVI(F2$) : LIMV=CVI(F3$)
140 CLOSE #1
150 DIM X(LIMC),Y(LIMC)
160 REM Coordinaten ...
170 OPEN "R",#1,"HDBCOORD.BIN",8
180 FIELD #1, 4 AS F1$, 4 AS F2$
190 FOR I=1 TO LIMC : GET #1,I
200 X(I)=CVS(F1$) : Y(I)=CVS(F2$)
210 REM PRINT X(I),Y(I)
220 NEXT I : CLOSE #1
230 REM
240 REM Plot mesh ...
250 SCREEN 2 : DIM NR(4)
260 OPEN "R",#1,"HDBTOPOL.BIN",16
270 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$
280 FOR I=1 TO LIMT : GET #1,I
290 NR(1)=CVI(F1$) : NR(2)=CVI(F2$)
295 NR(3)=CVI(F3$) : NR(4)=CVI(F4$)
300 FOR M=1 TO 4 : L=(M MOD 4)+1
310 X1=40*X(NR(M)) : Y1=190-20*Y(NR(M))
320 X2=40*X(NR(L)) : Y2=190-20*Y(NR(L))
330 LINE (X1,Y1)-(X2,Y2)
340 NEXT M : NEXT I : CLOSE #1
350 REM
360 REM Vectorplot snelheden ...
370 OPEN "R",#1,"HDBUNKNO.BIN",8
380 OPEN "R",#2,"HDBVELOC.BIN",8
390 FIELD #1, 4 AS F1$, 4 AS F2$
400 FIELD #2, 4 AS FU$, 4 AS FV$
410 FOR I=1 TO LIMV : GET #1,I : GET #2,I
420 N1=CVI(F1$) : N2=CVI(F2$) : UU=CVS(FU$) : VV=CVS(FV$)
430 X1=40*.5*(X(N1)+X(N2)) : Y1=190-20*.5*(Y(N1)+Y(N2))
440 X2=X1+40*UU : Y2=Y1-20*VV : LINE (X1,Y1)-(X2,Y2)
450 NEXT I : CLOSE #1 : CLOSE #2
460 END

10 REM Probleem van Labrujere
20 REM ======================
50 REM Snelheidsprofiel cylinder
60 REM -------------------------
80 REM Algemeen:
90 OPEN "R",#1,"HDBGENER.BIN",20
100 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$, 4 AS F5$
110 GET #1,1 : NPH=CVI(F5$)
120 CLOSE #1
130 DIM X(NPH+1),Y(NPH+1)
140 REM Coordinaten:
150 OPEN "R",#1,"HDBCOORD.BIN",8
160 FIELD #1, 4 AS F1$, 4 AS F2$
170 FOR I=1 TO NPH+1 : GET #1,I
180 X(I)=CVS(F1$) : Y(I)=CVS(F2$)
190 NEXT I : CLOSE #1
200 REM
210 XO=10 : YO=150 : SX=200 : SY=50 : SCREEN 2
220 LINE (XO,YO-SY)-(XO+SX,YO-SY) : LINE (XO,YO)-(XO,YO-4*SY)
230 REM
240 REM Plot snelheden:
250 OPEN "R",#1,"HDBVELOC.BIN",8
260 FIELD #1, 4 AS F1$, 4 AS F2$
270 FOR I=1 TO NPH : GET #1,I
280 U(I)=CVS(F1$) : V(I)=CVS(F2$)
290 NEXT I : CLOSE #1
300 REM
310 XN=XO+SX*X(1) : YN=YO-SY
320 REM Horizontale component:
330 FOR I=1 TO NPH-1
340 XM=XN : XN=XO+SX*.5*(X(I)+X(I+1))
350 YM=YN : YN=YO-SY*(1+U(I))
360 LINE(XM,YM)-(XN,YN) : NEXT I
370 LINE(XN,YN)-(XO+SX*X(NPH),YO-SY*(1+2))
380 XN=XO+SX*X(1) : YN=YO-SY
390 REM Verticale component:
400 XN=XO+SX*X(1) : YN=YO-SY
410 FOR I=1 TO NPH-1
420 XM=XN : XN=XO+SX*.5*(X(I)+X(I+1))
430 YM=YN : YN=YO-SY*(1+V(I))
440 LINE(XM,YM)-(XN,YN) : NEXT I
450 LINE(XN,YN)-(XO+SX*X(NPH),YO-SY)
460 REM
470 REM Analytische oplossingen:
480 LIMIT=25
490 XN=0 : YN=YO-SY
500 FOR I=1 TO LIMIT
510 PHI=2*ATN(1)/(LIMIT-1)*(I-1)
520 UU=1-COS(2*PHI)
530 XM=XN : XN=XO+COS(PHI)*SX
540 YM=YN : YN=YO-SY*(1+UU)
550 LINE(XM,YM)-(XN,YN) : NEXT I
560 XN=0 : YN=YO-SY
570 FOR I=1 TO LIMIT
580 PHI=2*ATN(1)/(LIMIT-1)*(I-1)
590 VV=-SIN(2*PHI)
600 XM=XN : XN=XO+COS(PHI)*SX
610 YM=YN : YN=YO-SY*(1+VV)
620 LINE(XM,YM)-(XN,YN) : NEXT I
630 END

10 REM Probleem van Labrujere
20 REM ======================
50 REM Stroomfunktie berekenen
60 REM -----------------------
90 OPEN "R",#1,"HDBGENER.BIN",20 : REM Algemeen
100 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$, 4 AS F5$
110 GET #1,1 : CLOSE #1
120 LIMC=CVI(F1$) : LIMV=CVI(F3$) : NPH=CVI(F5$)
130 NN=LIMC : NB1=NPH+1
140 REM
150 DIM X(LIMC),Y(LIMC),S(NN,NB1),B(NN)
160 REM
170 REM Coordinaten ...
180 OPEN "R",#1,"HDBCOORD.BIN",8
190 FIELD #1, 4 AS F1$, 4 AS F2$
200 FOR I=1 TO LIMC : GET #1,I
210 X(I)=CVS(F1$) : Y(I)=CVS(F2$)
220 NEXT I : CLOSE #1
230 REM
240 PRINT "Assemblage procedure voor";LIMV;"elementen ..."
250 REM
260 DIM NR(2),E(2,2),R(2)
270 FOR I=1 TO NN : FOR J=1 TO NB1 : S(I,J)=0
275 NEXT J : B(I)=0 : NEXT I
280 REM
290 OPEN "R",#1,"HDBUNKNO.BIN",8
295 FIELD #1, 4 AS F1$, 4 AS F2$
300 OPEN "R",#2,"HDBVELOC.BIN",8
305 FIELD #2, 4 AS FU$, 4 AS FV$
310 FOR F=1 TO LIMV : GET #1,F : GET #2,F : PRINT F;
320 NR(1)=CVI(F1$) : NR(2)=CVI(F2$)
325 U=CVS(FU$) : V=CVS(FV$)
330 REM
340 GOSUB 600 : REM Element matrix ...
350 REM
360 FOR I=1 TO 2 : IG=NR(I)
370 FOR J=I TO 2 : JG=NR(J)
380 II=IG : IF JG<IG THEN II=JG
390 JJ=ABS(IG-JG)+1
400 S(II,JJ)=S(II,JJ)+E(I,J)
410 NEXT J
420 B(IG)=B(IG)+R(I)
430 NEXT I
440 NEXT F : CLOSE #1 : PRINT
450 REM
460 REM Randvoorwaarde ...
470 S(NPH,1)=S(NPH,1)+10
480 REM
490 GOSUB 740 : REM Los vergelijkingen op ...
500 REM
510 OPEN "R",#1,"HDBPSI.BIN",4 : FIELD #1, 4 AS F1$
520 FOR K=1 TO LIMC
530 LSET F1$=MKS$(B(K)) : REM LPRINT B(K)
540 PUT #1,K : NEXT K : CLOSE #1
550 END
560 REM ------------------------
570 REM Eindige Elementen Matrix
580 REM ------------------------
590 REM Matrix schoonmaken:
600 FOR I=1 TO 2 : R(I)=0 : FOR J=1 TO 2 : E(I,J)=0
605 NEXT J : NEXT I
610 REM Differentie-schema:
620 X21=X(NR(2))-X(NR(1)) : Y21=Y(NR(2))-Y(NR(1))
630 A(1)=-1 : A(2)=+1 : W=Y21*U-X21*V
650 REM Kleinste kwadraten element:
660 FOR I=1 TO 2 : FOR J=I TO 2
670 E(I,J)=E(I,J)+A(I)*A(J)
680 NEXT J : R(I)=R(I)+A(I)*W : NEXT I
690 RETURN
700 REM  --------------------------------------------
710 REM  Oplossen van stelsel lineaire vergelijkingen
720 REM  Positief definiet symmetrische band matrices
730 REM  --------------------------------------------
740 IF NN<2 THEN PRINT "GAUSS: NN=";NN : END
750 IF NB1<2 THEN PRINT "GAUSS: NB1=";NB1 : END
760 REM
770 PRINT "Vegen voor";NN;"vergelijkingen ..."
780 II=NN-1 : KK=NB1-1
790 FOR I=1 TO II : DIAG=S(I,1) : PRINT I;
800 IF DIAG=0 THEN PRINT "Nul op diagonaal" : END
810 IF I>NN-NB1 THEN KK=NN-I
820 FOR K=1 TO KK
830 PIVOT=S(I,K+1)/DIAG : IF PIVOT=0 THEN GOTO 860
840 B(I+K)=B(I+K)-PIVOT*B(I) : JJ=KK-K+1
850 FOR J=1 TO JJ : S(I+K,J)=S(I+K,J)-PIVOT*S(I,K+J) : NEXT J
860 NEXT K : NEXT I : PRINT
870 PRINT "Terugsubstitutie ..."
880 KK=NB1-1 : FOR N=2 TO NN : I=NN-N+2 : PRINT I;
890 B(I)=B(I)/S(I,1) : IF I<NB1 THEN KK=I-1
900 FOR K=1 TO KK : B(I-K)=B(I-K)-B(I)*S(I-K,K+1) : NEXT K
910 NEXT N : B(1)=B(1)/S(1,1) : PRINT : RETURN
920 REM

10 REM Probleem van Labrujere
20 REM ======================
30 REM Contour plot van de stroom-funktie (stroomlijnen)
40 REM -------------------------------------------------
45 PRINT "Gegevens inlezen ..."
50 OPEN "R",#1,"HDBGENER.BIN",20
60 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$, 4 AS F5$
70 GET #1,1 : CLOSE #1 : LIMC=CVI(F1$) : LIMT=CVI(F2$)
90 DIM X(LIMC),Y(LIMC),B(LIMC)
110 REM
190 REM Coordinaten ...
200 OPEN "R",#1,"HDBCOORD.BIN",8
205 FIELD #1, 4 AS F1$, 4 AS F2$
210 FOR I=1 TO LIMC
220 GET #1,I : X(I)=CVS(F1$) : Y(I)=CVS(F2$)
230 NEXT I : CLOSE #1
240 REM
250 REM Oplossing ...
260 OPEN "R",#1,"HDBPSI.BIN",4 : FIELD #1, 4 AS F1$
270 FOR I=1 TO LIMC : GET #1,I : B(I)=CVS(F1$) : NEXT I
280 CLOSE #1
290 BOT=B(1) : TOP=B(1) : FOR I=2 TO LIMC
300 IF TOP<B(I) THEN TOP=B(I)
310 IF BOT>B(I) THEN BOT=B(I)
320 NEXT I
325 PRINT "Voorbereiding ..."
330 DIM ND(3,4)
332 DATA 1,2,3 , 1,2,4 , 4,3,1 , 4,3,2
334 FOR J=1 TO 4 : FOR I=1 TO 3 : READ ND(I,J)
335 NEXT I : NEXT J
336 DIM NR(4),XQ(4),YQ(4),ZQ(4),XD(3),YD(3),ZD(3),FD(3)
340 DIM NP(3,6),MUST(27),MOVE(27),R(6),NUM(3),NO(3)
350 GOSUB 1000
355 REM Doorloop mesh ...
360 SCREEN 2 : MOST=32 : N=MOST : GOSUB 920
370 OPEN "R",#1,"HDBTOPOL.BIN",16
380 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$
390 FOR F=1 TO LIMT : GET #1,F
400 NR(1)=CVI(F1$) : NR(2)=CVI(F2$)
405 NR(3)=CVI(F3$) : NR(4)=CVI(F4$)
410 FOR K=1 TO 4 : L=NR(K)
420 XQ(K)=X(L) : YQ(K)=Y(L) : ZQ(K)=(B(L)-BOT)/(TOP-BOT)
430 NEXT K
435 REM Vierhoeken:
440 QMIN=1 : QMAX=0
450 FOR K=1 TO 4 : REM Extrema:
460 IF ZQ(K)<QMIN THEN QMIN=ZQ(K)
470 IF ZQ(K)>QMAX THEN QMAX=ZQ(K)
480 NEXT K : REM Niets overbodigs:
490 Q=QMIN : GOSUB 800 : NQ1=N1
500 Q=QMAX : GOSUB 800 : NQ2=N2
505 REM
510 REM Doorloop driehoeken:
520 FOR D=1 TO 4 : FOR N=1 TO 3
530 K=ND(N,D) : XD(N)=XQ(K) : YD(N)=YQ(K) : ZD(N)=ZQ(K)
540 NEXT N
550 REM Doorloop niveaus:
560 FOR N=NQ1 TO NQ2 : ZERO=N/MOST
570 REM Elementaire contouren:
580 GOSUB 2000 : NEXT N : NEXT D
590 NEXT F : CLOSE #1 : END
790 REM
800 REM Bepaal onder & boven-grenzen
810 REM ----------------------------
830 N1=0 : N2=MOST : M=MOST
840 FOR K=1 TO LOG2+1
850 M=M\2 : N1=N1+M : N2=N2-M
860 IF Q<=N1/MOST THEN N1=N1-M
870 IF Q>=N2/MOST THEN N2=N2+M
880 NEXT K : RETURN
890 REM
900 REM Bepaal 2-log van getal
910 REM ----------------------
920 IF N<=0 THEN PRINT "LOG2: fout argument" : END
930 FOR K=1 TO 64
940 N=N\2 : LOG2=K-1 : IF N=0 THEN RETURN
950 NEXT K : PRINT "LOG2: draait rondjes" : END
1000 REM
1010 REM Initialiseer NP(3,6),MUST(27),MOVE(27)
1015 REM --------------------------------------
1020 REM Aantal elementen = 3
1030 REM Aantal eigenschappen = 3
1050 REM Maak alle permutaties van (0,1,2):
1060 FOR L=1 TO 6 : N=3 : GOSUB 1330
1070 FOR J=1 TO 3 : NP(J,L)=R(J) : NEXT J
1075 NEXT L
1080 REM Equivalentie-klassen:
1090 N=0 : FOR L=1 TO 27 : IF MUST(L)<>0 THEN 1220
1100 N=N+1 : REM Decimalen (in 3-tallig stelsel):
1110 M=L-1 : FOR K=1 TO 3
1115 NUM(3-K+1)=M MOD 3 : M=M\3 : NEXT K
1120 FOR KEER=1 TO 2
1125 FOR I=1 TO 6 : REM Inverse permutatie:
1130 FOR J=1 TO 3 : NO(NP(J,I))=J : NEXT J
1140 M=0 : FOR J=1 TO 3 : M=M*3+NUM(NO(J))
1145 NEXT J : M=M+1
1150 REM Toekennen van equivalente elementen:
1160 IF MUST(M)=0 THEN MUST(M)=N : REM Klasse
1170 IF MOVE(M)=0 THEN MOVE(M)=I : REM Permutatie
1180 NEXT I : REM Alle cijfer-permutaties behandeld.
1200 FOR K=1 TO 3 : NUM(K)=(3-NUM(K)) MOD 3 : NEXT K
1205 M=0 : FOR J=1 TO 3 : M=M*3+NUM(NO(J))
1207 NEXT J : M=M+1 : IF MUST(M)<>0 THEN 1220
1210 NEXT KEER : REM < en > de tweede KEER verwisseld
1220 NEXT L : RETURN
1260 REM
1270 REM Permutaties van array R(N)
1280 REM --------------------------
1290 REM L = volgnummer permutatie (invoer)
1300 REM N = aantal elementen van: (invoer)
1310 REM R = rij die de permutatie bevat (uitvoer)
1320 REM
1330 H=1 : FOR K=1 TO N : R(K)=K : H=H*K : NEXT K
1340 REM Welk item:
1350 FOR K=1 TO N : P=N-K+1 : H=H\P
1360 Q=(((L-1)\H) MOD P)+1 : IF Q=1 THEN 1400
1370 REM Roteer data:
1380 S=R(K+Q-1) : FOR M=K+Q-1 TO K+1 STEP -1
1390 R(M)=R(M-1) : NEXT M : R(K)=S
1400 NEXT K : RETURN
2000 REM
2010 REM Contourlijnen in een driehoek
2020 REM -----------------------------
2030 REM Invoer : XD(3),YD(3),ZD(3),ZERO
2050 FOR K=1 TO 3 : FD(K)=ZD(K)-ZERO : NEXT K
2060 REM
2065 REM Coderen van elk geval:
2070 MORE=0 : FOR K=1 TO 3
2080 IF FD(K)=0 THEN MORE=MORE*3
2090 IF FD(K)<0 THEN MORE=MORE*3+1
2100 IF FD(K)>0 THEN MORE=MORE*3+2
2110 NEXT K : LABEL=MUST(MORE+1)
2115 IF LABEL=1 OR LABEL=3 OR LABEL=5 THEN RETURN
2120 REM
2130 REM Permutatie van (1,2,3):
2140 KEUS=MOVE(MORE+1)
2145 NP1=NP(1,KEUS) : NP2=NP(2,KEUS) : NP3=NP(3,KEUS)
2150 X1=XD(NP1) : X2=XD(NP2) : X3=XD(NP3)
2160 Y1=YD(NP1) : Y2=YD(NP2) : Y3=YD(NP3)
2170 F1=FD(NP1) : F2=FD(NP2) : F3=FD(NP3)
2180 REM
2195 REM Aktie:     1    2    3    4    5    6
2200 ON LABEL GOTO 2210,2230,2210,2270,2210,2330
2210 REM Niets doen: === , =<< , <<<
2220 RETURN
2230 REM Lijn door EEN zijde: ==<
2240 LINE (50*X1,190-20*Y1)-(50*X2,190-20*Y2)
2260 RETURN
2270 REM Lijn door EEN hoekpunt: =<>
2275 ONDER=1/(F2-F3)
2280 XC=(F2*X3-F3*X2)*ONDER
2290 YC=(F2*Y3-F3*Y2)*ONDER
2300 LINE (50*X1,190-20*Y1)-(50*XC,190-20*YC)
2320 RETURN
2330 REM Lijn door twee zijden: <<>
2335 ONDER=1/(F1-F3)
2340 XB=(F1*X3-F3*X1)*ONDER
2350 YB=(F1*Y3-F3*Y1)*ONDER
2355 ONDER=1/(F2-F3)
2360 XC=(F2*X3-F3*X2)*ONDER
2370 YC=(F2*Y3-F3*Y2)*ONDER
2380 LINE (50*XB,190-20*YB)-(50*XC,190-20*YC)
2400 RETURN