BASIC NLR

10 REM Probleem van Labrujere
20 REM ======================
50 REM Mesh Generator
60 REM --------------
140 INPUT "Aantal rasterlijnen in radiale richting (7) : ",NDR
150 INPUT "Aantal rasterlijnen in omtreks-richting (5) : ",NPH
160 IF NDR<=0 THEN NDR=7
170 IF NPH<=0 THEN NPH=5
180 REM
190 REM ... Initialisatie ...
200 LIMC=NDR*NPH : REM aantal rasterpunten
205 LIMT=(NDR-1)*(NPH-1) : REM aantal elementen
210 LIMB=2*NDR+2*NPH : REM aantal randpunten
220 REM
230 PRINT "Algemene gegevens ..."
240 OPEN "R",#1,"NLRGENER.BIN",16
250 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$
260 LSET F1$=MKI$(LIMC) : LSET F2$=MKI$(LIMT)
270 LSET F3$=MKI$(LIMB) : LSET F4$=MKI$(NPH)
280 PUT #1,1 : CLOSE #1
290 REM
300 PRINT "Coordinaten ..."
310 OPEN "R",#1,"NLRCOORD.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 Buitengrens radiaal
350 FOR I=1 TO NDR
360 R=(1+P)^(I-1) : REM Geeft 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,"NLRTOPOL.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 "Randvoorwaarden ..."
580 OPEN "R",#1,"NLRBOUND.BIN",8
590 FIELD #1, 4 AS F1$, 4 AS F2$
600 FOR I=1 TO NDR : REM Instroming ...
610 NOB1=(I-1)*NPH+1 : NOB2=1
620 REM PRINT NOB1,NOB2
630 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
640 PUT #1,I : NEXT I : L=NDR
650 FOR I=1 TO NPH : REM Cylinder ...
660 NOB1=I : NOB2=2
670 REM PRINT NOB1,NOB2
680 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
690 PUT #1,I+L : NEXT I : L=L+NPH
700 FOR I=1 TO NDR : REM Symmetrie ...
710 NOB1=NPH*I : NOB2=1
720 REM PRINT NOB1,NOB2
730 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
740 PUT #1,I+L : NEXT I : L=L+NDR
750 FOR I=1 TO NPH : REM Uitstroming ...
760 NOB1=NPH*(NDR-1)+I : NOB2=3
770 REM PRINT NOB1,NOB2
780 LSET F1$=MKI$(NOB1) : LSET F2$=MKI$(NOB2)
790 PUT #1,I+L : NEXT I : CLOSE #1
800 END

10 REM Probleem van Labrujere:
20 REM ======================
30 REM Samenstellen en oplossen vergelijkingen
40 REM =======================================
50 REM Gegevens inlezen
60 REM ----------------
100 OPEN "R",#1,"NLRGENER.BIN",16
110 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$
120 GET #1,1 : LIMC=CVI(F1$) : LIMT=CVI(F2$) 
125 LIMB=CVI(F3$) : NPH=CVI(F4$)
130 CLOSE #1 : NN=2*LIMC : NB1=2*NPH+4
140 REM
150 DIM X(LIMC),Y(LIMC),S(NN,NB1),B(NN)
160 REM
170 REM Coordinaten:
180 OPEN "R",#1,"NLRCOORD.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 REM Grote matrix en vector schoonmaken:
250 FOR I=1 TO NN : FOR J=1 TO NB1 : S(I,J)=0 : NEXT J
255 B(I)=0 : NEXT I
260 REM
270 DIM NO(4),NR(8),E(8,8),R(2),A(6),NG(3,2)
280 REM
290 OPEN "R",#1,"NLRTOPOL.BIN",16
300 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$ 
305 GET #1,1
310 NO(1)=CVI(F1$) : NO(2)=CVI(F2$) 
315 NO(3)=CVI(F3$) : NO(4)=CVI(F4$)
320 PRINT "Repeteerbaar element ..." : GOSUB 830
330 REM
340 PRINT "Assemblage procedure voor";LIMT;"elementen ..."
360 FOR F=1 TO LIMT : GET #1,F : PRINT F;
370 NO(1)=CVI(F1$) : NO(2)=CVI(F2$)
375 NO(3)=CVI(F3$) : NO(4)=CVI(F4$)
380 REM Nummering snelheden:
390 FOR I=1 TO 4 : FOR J=1 TO 2 : K=2*(I-1)+J
400 NR(K)=2*(NO(I)-1)+J : NEXT J : NEXT I
410 REM
420 REM Bulk elementen intellen:
430 FOR I=1 TO 8 : IG=NR(I)
440 FOR J=I TO 8 : JG=NR(J)
450 II=IG : IF JG<IG THEN II=JG
460 JJ=ABS(IG-JG)+1
470 S(II,JJ)=S(II,JJ)+E(I,J)
480 NEXT J : NEXT I
490 NEXT F : CLOSE #1 : PRINT
500 REM
510 PRINT "Randvoorwaarden ..."
520 OPEN "R",#1,"NLRBOUND.BIN",8
530 FIELD #1, 4 AS F1$, 4 AS F2$
540 FOR F=1 TO LIMB : GET #1,F
550 MU=CVI(F1$) : KIND=CVI(F2$)
560 GOSUB 1780
570 REM Nummering snelheden:
580 NR(1)=2*MU-1 : NR(2)=2*MU
590 REM Elementen intellen:
600 FOR I=1 TO 2 : IG=NR(I)
610 FOR J=I TO 2 : JG=NR(J)
620 II=IG : IF JG<IG THEN II=JG
630 JJ=ABS(IG-JG)+1
640 S(II,JJ)=S(II,JJ)+E(I,J)
650 NEXT J
660 B(IG)=B(IG)+R(I)
670 NEXT I
680 NEXT F : CLOSE #1 : REM Einde assemblage
690 REM
700 GOSUB 1980 : REM Vergelijkingen oplossen
710 REM
720 OPEN "R",#1,"NLRVELOC.BIN",8
730 FIELD #1, 4 AS F1$, 4 AS F2$
740 FOR K=1 TO LIMC
750 LSET F1$=MKS$(B(2*K-1)) : LSET F2$=MKS$(B(2*K))
760 REM LPRINT B(2*K-1),B(2*K)
770 PUT #1,K : NEXT K : CLOSE #1
780 END
790 REM ==========================
800 REM Eindige elementen matrices
810 REM ==========================
820 REM Maak element schoon:
830 FOR I=1 TO 8 : FOR J=1 TO 8 : E(I,J)=0 : NEXT J : NEXT I
910 NG(1,1)=1 : NG(2,1)=2 : NG(3,1)=4
920 NG(1,2)=2 : NG(2,2)=3 : NG(3,2)=4
930 FOR IE=1 TO 2
940 REM ... Geometrie ...
950 NO1=NO(NG(1,IE)) : NO2=NO(NG(2,IE)) : NO3=NO(NG(3,IE))
960 X21=X(NO2)-X(NO1) : Y21=Y(NO2)-Y(NO1)
970 X32=X(NO3)-X(NO2) : Y32=Y(NO3)-Y(NO2)
980 X13=X(NO1)-X(NO3) : Y13=Y(NO1)-Y(NO3)
990 GOSUB 1230 : NEXT IE : RETURN
1190 REM ------------------------------
1200 REM Eindige Differentie benadering
1210 REM ------------------------------
1220 REM ... Onsamendrukbaar ...
1230 A(1)=-Y32 : A(2)=+X32 : A(3)=-Y13 
1235 A(4)=+X13 : A(5)=-Y21 : A(6)=+X21
1240 GOSUB 1580 : REM Voeg toe aan element ...
1250 REM ... Rotatievrij ...
1260 A(1)=+X32 : A(2)=+Y32 : A(3)=+X13 
1265 A(4)=+Y13 : A(5)=+X21 : A(6)=+Y21
1270 GOSUB 1580 : REM Voeg toe aan element ...
1280 RETURN
1560 REM --------------------------
1570 REM Kleinste Kwadraten element
1580 REM --------------------------
1590 REM SpoorWegen:
1600 SUM=0 : FOR I=1 TO 6 : SUM=SUM+A(I)*A(I)
1605 NEXT I : SUM=1/SUM
1610 REM Eindig element:
1620 FOR I=1 TO 3 : FOR K=1 TO 2
1630 IH=2*(I-1)+K : IG=2*(NG(I,IE)-1)+K
1640 FOR J=1 TO 3 : FOR L=1 TO 2
1650 JH=2*(J-1)+L : JG=2*(NG(J,IE)-1)+L
1660 E(IG,JG)=E(IG,JG)+A(IH)*A(JH)*SUM
1670 NEXT L : NEXT J : NEXT K : NEXT I
1680 RETURN
1760 REM ===============
1770 REM Randvoorwaarden
1780 REM ===============
1790 FOR I=1 TO 2 : R(I)=0 : FOR J=1 TO 2
1795 E(I,J)=0 : NEXT J : NEXT I
1800 REM
1810 ON KIND GOTO 1840,1870,1920
1820 REM
1830 REM ... V = 0 ...
1840 E(2,2)=1 : RETURN
1850 REM
1860 REM ... X.U + Y.V = 0 ...
1870 X1=X(MU) : Y1=Y(MU) : SUM=X1*X1+Y1*Y1
1880 E(1,1)=X1*X1/SUM : E(2,2)=Y1*Y1/SUM
1890 E(1,2)=X1*Y1/SUM : E(2,1)=X1*Y1/SUM
1900 RETURN
1910 REM ... U = 1 ...
1920 E(1,1)=1 : R(1)=1 : RETURN
1930 REM
1940 REM ==========================================
1950 REM Oplossen van linear stelsel vergelijkingen
1960 REM Positief definiete symmetrische bandmatrix
1970 REM ==========================================
1980 IF NN<2 THEN PRINT "GAUSS: NN=";NN : END
1990 IF NB1<2 THEN PRINT "GAUSS: NB1=";NB1 : END
2000 REM
2010 PRINT "Vegen voor";NN;"vergelijkingen ..."
2020 II=NN-1 : KK=NB1-1
2030 FOR I=1 TO II : DIAG=S(I,1) : PRINT I;
2040 IF DIAG=0 THEN PRINT "Nul op diagonaal" : END
2050 IF I>NN-NB1 THEN KK=NN-I
2060 FOR K=1 TO KK
2070 PIVOT=S(I,K+1)/DIAG : IF PIVOT=0 THEN GOTO 2100
2080 B(I+K)=B(I+K)-PIVOT*B(I) : JJ=KK-K+1
2090 FOR J=1 TO JJ : S(I+K,J)=S(I+K,J)-PIVOT*S(I,K+J) : NEXT J
2100 NEXT K : NEXT I
2110 PRINT "Terugsubstitutie ..."
2120 KK=NB1-1 : FOR N=2 TO NN : I=NN-N+2 : PRINT I;
2130 B(I)=B(I)/S(I,1) : IF I<NB1 THEN KK=I-1
2140 FOR K=1 TO KK : B(I-K)=B(I-K)-B(I)*S(I-K,K+1) : NEXT K
2150 NEXT N : B(1)=B(1)/S(1,1) : RETURN
2160 REM

10 REM Probleem van Labrujere
20 REM ======================
50 REM MESH PLOT
60 REM ---------
100 REM Algemeen ...
110 OPEN "R",#1,"NLRGENER.BIN",16
120 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$
130 GET #1,1 : LIMC=CVI(F1$) : LIMT=CVI(F2$)
140 CLOSE #1
150 DIM X(LIMC),Y(LIMC)
160 REM Coordinaten ...
170 OPEN "R",#1,"NLRCOORD.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,"NLRTOPOL.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 ...
380 OPEN "R",#2,"NLRVELOC.BIN",8
400 FIELD #2, 4 AS FU$, 4 AS FV$
410 FOR I=1 TO LIMC : GET #2,I
420 UU=CVS(FU$) : VV=CVS(FV$)
430 X1=40*X(I) : Y1=190-20*Y(I)
440 X2=X1+40*UU : Y2=Y1-20*VV : LINE (X1,Y1)-(X2,Y2)
450 NEXT I
460 END

10 REM Probleem van Labrujere
20 REM ======================
50 REM CYLINDER SNELHEIDS PROFIEL
60 REM --------------------------
80 REM Algemeen ...
90 OPEN "R",#1,"NLRGENER.BIN",16
100 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$, 4 AS F4$
110 GET #1,1 : NPH=CVI(F4$)
120 CLOSE #1
130 DIM X(NPH),Y(NPH),U(NPH),V(NPH)
140 REM Coordinaten ...
150 OPEN "R",#1,"NLRCOORD.BIN",8
160 FIELD #1, 4 AS F1$, 4 AS F2$
170 FOR I=1 TO NPH : 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,"NLRVELOC.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=2 TO NPH
340 XM=XN : XN=XO+SX*X(I)
350 YM=YN : YN=YO-SY*(1+U(I))
360 LINE(XM,YM)-(XN,YN) : NEXT I
390 REM Verticale component ...
400 XN=XO+SX*X(1) : YN=YO-SY
410 FOR I=2 TO NPH
420 XM=XN : XN=XO+SX*X(I)
430 YM=YN : YN=YO-SY*(1+V(I))
440 LINE(XM,YM)-(XN,YN) : NEXT I
460 REM
470 REM Analytische oplossing ...
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