' Terbeli pontokra "legkozelebbi" sikot illeszt. (Regresszioszamitas.)
' Leiras a program vegen.

RANDOMIZE TIMER

mintav = 1000000

SCREEN 13

'szinek
PALETTE 0, 0
PALETTE 250, 65536 * 63 + 256 * 63 + 63
PALETTE 251, 63
PALETTE 252, 65536 * 63 + 63
PALETTE 253, 65536 * 11 + 256 * 11 + 11

FOR i = 1 TO 240
IF i < 60 THEN a = 60: b = i: c = 0
IF i < 120 AND i >= 60 THEN b = 60: a = 120 - i: c = 0
IF i >= 120 AND i <= 180 THEN a = 0: b = 60: c = i - 120
IF i >= 180 AND i < 210 THEN a = 0: b = 240 - i: c = b
IF i >= 210 THEN a = i - 210: b = 30: c = 30
PALETTE i, 65536 * a + 256 * b + c
NEXT i

n = 11
DIM vx(n)
DIM vy(n)
DIM vz(n)
pi = 3.14159265#

elore:
CLS

FOR i = 1 TO 240
LINE (i + 10, 196)-(i + 10, 200), i
NEXT i

'Veletlenszeru pontokat general.

sx = 0: sy = 0: sz = 0
sx2 = 0: sy2 = 0
sxy = 0: sxz = 0: syz = 0

FOR i = 1 TO n
'tetaa = (5 * RND + 5) * pi / 21
'tetaa = (RND - .5) * pi
tetaa = (RND - .5) * pi / 2
'tetaa = pi / 7
fii = i * 2 * pi / n
vx(i) = 20 * COS(tetaa) * COS(fii)
vy(i) = 20 * COS(tetaa) * SIN(fii)
vz(i) = 20 * SIN(tetaa)
NEXT i

'GOTO nemtokveletlen
' Ha ez a "goto" mukodik, akkor az aposztroffal nem vedett kis tetakra
' (+-pi/6) a kozelites eleg jo lesz. Ha meg nem mukodik (vagyis aposztroffal
' van vedve), akkor rossz lesz a kozelites, de sokkal szebb az abra.

FOR i = 1 TO n
vx(i) = RND * 50 - 25
vy(i) = RND * 50 - 25
vz(i) = RND * 50 - 25
NEXT i

nemtokveletlen:

FOR k2 = 0 TO 150 STEP 25
LINE (1, k2 + 5)-(5, k2 + 5), 250
LINE (306, k2 + 5)-(310, k2 + 5), 250
NEXT k2

FOR k1 = 0 TO 300 STEP 25
LINE (k1 + 5, 1)-(k1 + 5, 5), 250
LINE (k1 + 5, 156)-(k1 + 5, 160), 250
NEXT k1

'Atlagot szamol.

svx = 0: svy = 0: svz = 0

FOR i = 1 TO n
svx = svx + vx(i)
svy = svy + vy(i)
svz = svz + vz(i)
NEXT i

svx = svx / n
svy = svy / n
svz = svz / n

' A konnyebb megoldhatosag erdekeben csak olyan pontokat vizsgal, amelyek
' sulypontja az origo. Ennek erdekeben modositja a kezdeti, veletlenszeruen
' generalt pontokat.

FOR i = 1 TO n
vx(i) = vx(i) - svx
vy(i) = vy(i) - svy
vz(i) = vz(i) - svz

'A 3. kozeliteshez kepezi a szukseges adatokat.
sx = sx + vx(i)
sy = sy + vy(i)
sz = sz + vz(i)
sx2 = sx2 + vx(i) ^ 2
sy2 = sy2 + vy(i) ^ 2
sxy = sxy + vx(i) * vy(i)
sxz = sxz + vx(i) * vz(i)
syz = syz + vy(i) * vz(i)
NEXT i

FOR i1 = -75 TO 75
teta = pi / 150 * i1
FOR i2 = 1 TO 300
fi = pi / 150 * i2

'Egysegvektort kepez; ezzel mint normalvektorral jellemzi az illeszteni
'kivant sikot.

rx = COS(teta) * COS(fi)
ry = COS(teta) * SIN(fi)
rz = SIN(teta)

tav = 0

FOR i = 1 TO n
'tav = tav + ABS(vx(i) * rx + vy(i) * ry + vz(i) * rz)
tav = tav + (vx(i) * rx + vy(i) * ry + vz(i) * rz) ^ 2 / 20
NEXT i

IF i2 > 150 THEN GOTO nemmin
IF tav < mintav THEN mintav = tav: mini1 = i1: mini2 = i2
nemmin:
PSET (i2 + 5, -i1 + 80), 1 + INT(tav - 240 * INT(tav / 240))
NEXT i2
NEXT i1

CIRCLE (mini2 + 5, -mini1 + 80), 2, 251

svx = 0: svy = 0: svz = 0

'1. kozelites: a vektorok osszege volna a regresszios sik normalvektora.
'Ez nem hasznalhato erdemben, miutan a vektorokat olyanokka konvertaltuk,
'hogy osszeguk nullvektor legyen.

'FOR i = 1 TO n
'svx = svx + vx(i)
'svy = svy + vy(i)
'svz = svz + vz(i)
'NEXT i

'hossz = SQR(svx ^ 2 + svy ^ 2 + svz ^ 2)
'IF hossz = 0 THEN hossz = .1
'teta1 = ATN((SQR(svx ^ 2 + svy ^ 2)) / hossz)
'IF svx = 0 THEN svx = .1
'fi1 = ATN(svy / svx)
'IF fi1 < 0 THEN fi1 = fi1 + pi
'CIRCLE (fi1 / pi * 150 + 5, -teta1 / pi * 150 + 80), 1, 250


'2. kozelites: a vektorok sorbarakasa utan parosaval vesszuk a vektori
'szorzatokat, es ezek osszege volna a regresszios sik normalvektora.

svx = 0: svy = 0: svz = 0
FOR i = 1 TO n - 1
svx = svx + vy(i) * vz(i + 1) - vz(i) * vy(i + 1)
svy = svy + vz(i) * vx(i + 1) - vx(i) * vz(i + 1)
svz = svz + vx(i) * vy(i + 1) - vy(i) * vx(i + 1)
NEXT i

svx = svx + vy(n) * vz(1) - vz(n) * vy(1)
svy = svy + vz(n) * vx(1) - vx(n) * vz(1)
svz = svz + vx(n) * vy(1) - vy(n) * vx(1)

' Polarszogeket szamol a minimalizalonak szant xyz adatokbol.

IF svx ^ 2 + svy ^ 2 < .001 AND svz > 0 THEN teta2 = pi / 2: fi2 = pi: SOUND (1000), .1: GOTO nemarctg
IF svx ^ 2 + svy ^ 2 < .001 AND svz <= 0 THEN teta2 = -pi / 2: fi2 = pi: SOUND (1000), .1: GOTO nemarctg
teta2 = ATN(svz / SQR(svx ^ 2 + svy ^ 2))

IF ABS(svx) < .001 AND svy >= 0 THEN fi2 = pi / 2: SOUND (1000), .1: GOTO nemarctg3
IF ABS(svx) < .001 AND svy < 0 THEN fi2 = -pi / 2: SOUND (1000), .1: GOTO nemarctg3
fi2 = ATN(svy / svx)

IF svx < 0 AND svy > 0 THEN fi2 = fi2 + pi
IF svx < 0 AND svy < 0 THEN teta2 = -teta2
IF svx > 0 AND svy < 0 THEN fi2 = fi2 + pi: teta2 = -teta2

nemarctg:

CIRCLE (fi2 / pi * 150 + 5, -teta2 / pi * 150 + 80), 1, 252

minteta = pi / 150 * mini1
minfi = pi / 150 * mini2

cosiv = COS(pi / 2 - minteta) * COS(pi / 2 - teta2) + SIN(pi / 2 - minteta) * SIN(pi / 2 - teta2) * COS(fi2 - minfi)
siniv = SQR(1 - cosiv ^ 2)
ivelteres = ABS(ATN(siniv / cosiv))
COLOR 252
LOCATE 21, 1: PRINT "2. tevedes: "; ivelteres

' 3. kozelites: a kincstari regresszios sik meredeksegeet szamolja ki, ebbol
' normalvektort kepez, ebbol meg tavolsagot szamol. Ellenorzi, hogy a suly-
' ponton atmeno sikok kozul is ez-e a tavolsagminimalizalo allasu sik. (Nem.)

deta = n * sx2 * sy2 + 2 * sxy * sx * sy - sy2 * sx ^ 2 - sx2 * sy ^ 2 - n * sxy ^ 2
detAa = n * sxz * sy2 + sxy * sy * sz + syz * sx * sy - sy2 * sx * sz - sxz * sy ^ 2 - n * sxy * syz
detAb = n * sx2 * syz + sxz * sx * sy + sxy * sx * sz - syz * sx ^ 2 - sx2 * sy * sz - n * sxy * sxz
detAc = sx2 * sy2 * sz + sxy * syz * sx + sxz * sxy * sy - sy2 * sxz * sx - sx2 * sy * syz - sz * sxy ^ 2

IF detAa ^ 2 + detAb ^ 2 < .001 AND deta < 0 THEN teta3 = pi / 2: fi3 = pi: SOUND (1000), .1: GOTO nemarctg3
IF detAa ^ 2 + detAb ^ 2 < .001 AND deta >= 0 THEN teta3 = -pi / 2: fi3 = pi: SOUND (1000), .1: GOTO nemarctg3
teta3 = ATN(-deta / SQR(detAa ^ 2 + detAb ^ 2))

IF ABS(detAa) < .001 AND detAb >= 0 THEN fi3 = pi / 2: SOUND (1000), .1: GOTO nemarctg3
IF ABS(detAa) < .001 AND detAb < 0 THEN fi3 = -pi / 2: SOUND (1000), .1: GOTO nemarctg3
fi3 = ATN(detAb / detAa)

IF detAa < 0 AND detAb > 0 THEN fi3 = fi3 + pi
IF detAa < 0 AND detAb < 0 THEN teta3 = -teta3
IF detAa > 0 AND detAb < 0 THEN fi3 = fi3 + pi: teta3 = -teta3

nemarctg3:
CIRCLE (fi3 / pi * 150 + 5, -teta3 / pi * 150 + 80), 1, 253

minteta = pi / 150 * mini1
minfi = pi / 150 * mini2

'Gombi costetel
cosiv = COS(pi / 2 - minteta) * COS(pi / 2 - teta3) + SIN(pi / 2 - minteta) * SIN(pi / 2 - teta3) * COS(fi3 - minfi)
siniv = SQR(1 - cosiv ^ 2)
ivelteres3 = ABS(ATN(siniv / cosiv))
COLOR 253
LOCATE 22, 1: PRINT "3. tevedes: "; ivelteres3

dd: d$ = INKEY$
IF d$ = "" THEN GOTO dd

'4. kozelites: a z koordinatakbeli kulonbsegek negyzetosszege lesz a norma.

'Ujrarajzolja a domborzati terkepet e masik norma szerint. Ez az, ami alapjan
'a harmadik kozelites is dolgozik, es ennek itteni ervenyesseget van hivatva
'megmutatni a mostani abrazolas.

'A kijelzes a jobb oldali felkepernyore tortenik, es nem azt a teruletet
'abrazolja , amelyet most felulir, hanem ugyanazt, ami a bal felkepernyon
'lathato. Igy lehet jol osszehasonlitani, hogy a ketfele norma szerint nem
'ugyanott vannak a dombok tetoi, sem a godrok fenekei (igaz, a maximumot
'senki nem keresi). Mutatja azt is, hogy a z-"tavolsag" a vegtelenbe szall,
'ha a sik fuggolegeshez tart.

minsdz2 = 1000000

LINE (155, 1)-(155, 155), 250
LINE (156, 5)-(305, 155), 0, BF
FOR i1 = -75 TO 75
teta = pi / 150 * i1
FOR i2 = 1 TO 150
fi = pi / 150 * i2

rx = COS(teta) * COS(fi)
ry = COS(teta) * SIN(fi)
rz = SIN(teta)
IF ABS(rz) < .01 THEN GOTO nemmered
meredx = -rx / rz
meredy = -ry / rz

sdz2 = 0

FOR i = 1 TO n
'masikz = meredx * vx(i) + meredy * vy(i) + detAc / detA
masikz = meredx * vx(i) + meredy * vy(i)
sdz2 = sdz2 + (vz(i) - masikz) ^ 2 / 20
NEXT i

IF sdz2 < minsdz2 THEN minsdz2 = sdz2: mini2 = i2: mini1 = i1
PSET (i2 + 155, -i1 + 80), 1 + INT(sdz2 - 240 * INT(sdz2 / 240))
'LOCATE 1, 1: PRINT sdz2

nemmered:
NEXT i2
NEXT i1

CIRCLE (mini2 + 155, -mini1 + 80), 2, 251
CIRCLE (fi3 / pi * 150 + 155, -teta3 / pi * 150 + 80), 1, 253

minteta = pi / 150 * mini1
minfi = pi / 150 * mini2

cosiv = COS(pi / 2 - minteta) * COS(pi / 2 - teta3) + SIN(pi / 2 - minteta) * SIN(pi / 2 - teta3) * COS(fi3 - minfi)
siniv = SQR(1 - cosiv ^ 2)
ivelteres = ABS(ATN(siniv / cosiv))
COLOR 253: LOCATE 23, 1: PRINT "4. tevedes: "; ivelteres

ee: e$ = INKEY$
IF e$ = "" THEN GOTO ee
IF e$ <> CHR$(13) THEN COLOR 250: END
mintav = 1000000
GOTO elore

' **********************************************************************

' Veletlenszeruen megadunk n db pontot a terben ugy, hogy sulypontjuk az
' origo legyen. Keressuk azt az origon atmeno sikot, amelytol merve a pontok
' tavolsagainak negyzetosszege minimalis. A keresest a sik normalvektoranak
' gombi koordinatak szerinti pasztazo mozgatasaval vegezzuk.

' Az abrazolas lenyege: a polarszogek (teta es fi) ciklusos vegigjarasaval
' minden lehetseges sikallasra (normalvektor-allasra) meghatarozza az illeto
' siktol valo tavolsagnegyzet-osszeget, s az igy kapott erteket a terkepi
' magassagi szinezeshez hasonlo modon kepernyore irja. Az abra ket fele
' csusztatva tukrozessel egymasba viheto, de a kettevagassal gyakran eppen
' a minimum- es maximumhelyek domborzati viszonyai csonkulnanak meg, ezert
' inkabb ketszer abrazoltam minden iranyt.

' A minimalizalo sikallast a program gepiesen megkeresi: ez a pontos ertek.

' Ezzel a gep osszehasonlitja a folvett pontok helyvektoraibol valamilyen
' kozelito eljaras vagy keplet segitsegevel kapott ("sejtett") minimalizalo
' normalvektort, vagyis az e modszer szerinti regresszios sikot.

' Az elso, durva kozelitesi modszer (a vektorok osszege) nem mukodik, ha ez
' nullvektor.

' A masodik kozelitesi modszer (a vektorok onkenyes sorbarendezesevel
' letrehozott szomszedokbol kepezett vektori szorzatok osszege) csak
' azokra a pont-elhelyezesekre mukodik turhetoen, amelyekben a pontok
' a gomb valamely fo"korehez (az egyenlitohoz) kozel esnek (+-45 fok).

' A harmadik modszer nem mas, mint a kincstari regresszioszamito keplet,
' amely azonban nagy szamitasigenyu, azaz lassu.

' A negyedik modszer megegyezik a harmadikkal, de a tavolsagnegyzetek
' helyett a z koordinatakban fennallo kulonbsegek negyzetosszege alapjan
' minositi a sikillesztest jonak vagy rossznak.

' A vizszintes tengelyen az oraszog (fi) megy 0-tol 2pi-ig, majd a kepernyo
' ketteosztasa utan 0-tol pi-ig mindket felkepernyon.
' A fuggoleges tengelyen lentrol folfele a foldrajzi szelesseg szoge (teta)
' megy -pi/2-tol pi/2-ig.

' A sarkpontokat ez az abrazolas foloslegesen szethuzza az also es felso
' vonalakka, de ez masodlagos kerdes. Fo, hogy font es lent a vizszintes
' elteresek nem szamitanak. A prg. a "tevedest" (ivbeli elterest) a gombi
' koszinusztetellel kiszamitja es kiirja. Ehhez tudni kell, hogy az abra
' aljan es tetejen levo "pontok" igen kozel is lehetnek egymashoz, hiszen az
' altaluk kepviselt normalvektoroknak csak az allasa szamit, az iranya nem.

' Ertekeles: Csakis a 4. kozelites megbizhato mindene esetben, mely a z
' szerinti "legkisebb negyzetek modszere" segitsegevel mukodik.
' Papiron hasztalan probaltam levezetni olyan kepletet, amely a sikra
' meroleges tavolsag negyzetere explicite megadja az illeszkedo sikallast.
' Akar Descartes-, akar gombi polarkoordinatakban probalkoztam, az
' ismeretlenek csunyan fuggtek a pontok koordinatainak osszegetol.
' Illetve a helyvektorok osszege jott ki optimalis normalvektorul, de
' alighanem hibas a levezetes.

' Tavlati cel: idotakarekos kozelito modszer valamely terbeli ponthalmaz
' erintosikjanak meghatarozasara.

' Meg tavolibb cel: negydimenzios fraktalok abrazolasa mozgo haromdimenzios
' animacio latszati kepekent. Ehhez az elobbi eredmeny a latszo feluleti
' fenyesseg meghatarozasa vegett szukseges (onarnyek, esetleg csillogas).

' 2007. IV.29. Nemeth Ferenc (nemo44@hotmail.com)