Programma in QBasic difettoso

di il
21 risposte

21 Risposte - Pagina 2

  • Re: Programma in QBasic difettoso

    Elimina la

    On Error Goto

    che hai inserito e che esegue solo una

    Resume Next

    e così non serve proprio a nulla.

    Ti verrà evidenziato l'errore che è un “Subscript out of range” ovvero i valori di N1 o N2 nella funzione non sono corretti.

  • Re: Programma in QBasic difettoso

    Qui https://visualstudio.microsoft.com/it/vs/community/ trovi visual studio community 

  • Re: Programma in QBasic difettoso

    La simulazione degli N-corpi e' uno di quei problemi che chiunque voglia definirsi “programmatore” dovrebbe affrontare almeno una volta nella vita ;-)

    Il tuo codice mi sembra inutilmente complicato. Forse varebbe la pena semplificarlo!

    L'algo per simulare il comportamento di N corpi e' molto semplice, non servono “sofisticate” conoscenze di programmazione, ma PURA programmazione procedurale (il tipo di programmazione piu' semplice) ;-)

    Il problema lo puoi rappresentare con SOLO 3 oggetti: sia N il numero di corpi

    1. un vettore di lunghezza N con le masse dei corpi
    2. una matrice Nx6 contenente: nelle prime 3 colonne il VETTORE POSIZIONE dei corpi, e nelle seconde 3 colonne il VETTORE VELOCITA all'istante t
    3. un'altra matrice Nx6 che conterra'  POSIZIONE e VELOCITA' all'istante (t+1)
    4. poi scambierai le due matrici. NON SERVE fare delle copie ma SOLO usarle al contrario. Il trucco sta' nel cosiderarer la parita' di t (t, un intero, pari o dispari)

    Poi ti serve un po' di trigonometrica per proiettare la posizione dei corpi sul piano DATO il punto di vista in 3 dimensioni. Queste sono semplici formulette di grafica 3D (SEMPLICI, NON banali ;-) )

    Il problema ha complessita' N^2, cioe' il numero di calcoli da fare cresce con il quadrato del numero di corpi. Con 100 corpi fai piu' o meno 10_000 operazioni “complesse”. Il problema e' che con 1000 corpi te ne servono 1000_000 e con 1000_0000 di corpi te ne servono 10^12 (parecchio tantini ;-))

    La parte divertente e' iniziare a lavorare con 1000/10_000 di corpi ;-) 
    I mille (e forse anche i diecimila) corpi si possono calcolare anche con un personal computer di ultima generazione (diciamo ultimi 5/10 anni) 

    MA 

    per gestirli senza spendere settimane servono algo “intelligenti” ;-).
    Qui NON E' un problema di codifica MA di algoritmo ;-)

    Ed inventare algo intelligenti non ha niente a che fare con la programmazione. 

    Un'altra cosa da fare, e' anche quella di cambiare linguaggio di programmazione ed usare il C (non serve il C++babasta il C puro) , in grado di sfruttare meglio l'hardware a disposizione (come si suol dire: grattare direttamente il ferro ;-) ) Comunque questo impatta solo sul numero massimo di corpi simula ile in tempi ragionevoli. 

    Con un po' di pazienza trovi parecchia documentazione al riguardo: cerca “N body simulation | problem | software”. 
    Ed anche algoritmi gia' pronti (per simulazioni “complesse”).

    Puoi iniziare da qui:

    https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation

    https://en.wikipedia.org/wiki/N-body_simulation

  • Re: Programma in QBasic difettoso

    Tra l'altro, analizzando il tuo codice (i 2 cicli for innestati) hai fatto un errore grossolano! 

    sia i l'i-mo corpo. per calcolare posizione e velocità di i all'istante t+1, stai usando posizione e velocità all'istante t+1 per i corpi j<i e all'istante t per i corpi j>i.

    “Questo calcolo non sa da fare né oggi né mai” come disse qualcuno qualche tempo fa a proposito di un matrimonio ;-) 

  • Re: Programma in QBasic difettoso

    Vi ringrazio tutti per i consigli che mi state dando, cercherò di farne tesoro.

  • Re: Programma in QBasic difettoso

    Grazie ai vostri consigli ho migliorato moltissimo il programma, perciò ho deciso di condividerne con voi la versione (credo) definitiva. Avevo scelto di non utilizzare gli algoritmi semplificatori che mi avete proposto perché volevo scoprire fino a che punto potessi tenere insieme precisione dei calcoli e fluidità di esecuzione, poi alla fine ho ridotto il numero massimo dei corpi a 10; mi piacerebbe sapere se “traducendolo” in un altro linguaggio riuscirei ad aumentare N senza rallentarlo troppo, o se la complessità del problema - che, se non ho capito male, cresce con N^2 - sia indipendente dal linguaggio di programmazione.
    P.S.: poiché mi avete spiegato che ho sbagliato sezione proponendo questo thread qui, potreste dirmi in quale sezione vanno esposti i problemi relativi ai programmi in QBasic, visto che ne ho un altro da postare? Grazie.

    $Debug
    ' NEWTON
    ' Simulazione del moto di molti corpi celesti
    ' dimensiona e definisci funzioni
    Dim G As Double, VC As Double, AL As Double
    Dim it As Double, pr As Double, dt As Double, pigreco As Double
    Dim nc As Integer, ncmax As Integer, Nt As Integer
    Dim N1 As Integer, N2 As Integer
    Dim MT As Double, RT As Double, Dsatg As Double
    Dim MS As Double, RS As Double, DTerra As Double
    Dim Mtot As Double, Mmax As Integer, Dmax As Integer
    Dim vmax1 As Double, vmax2 As Double
    Dim q0 As Double, dq As Single
    Dim fndist As Double, fndistmax As Double
    Dim Quota As Double, QApo As Double, QPer As Double
    Dim acc As Double
    Dim k As Double, L As Integer, C As Integer
    Dim lo As Integer, lv As Integer
    Dim pm As Integer, primo As Integer
    Dim I1 As Integer, I2 As Integer
    Dim diag1 As Integer, diag2 As Integer, diag3 As Integer
    Dim AK As Integer, AP As Single, AN As Single, DA As Single
    Dim scr As Integer, col As Integer
    Dim m1 As Integer, m2 As Integer, den As Single, Ecs As Double
    Dim Ep As Double, Epmax As Double, Epmin As Double, Epmed As Double
    Dim Ec As Double, Ecmax As Double, Ecmin As Double, Ecmed As Double
    Dim Etot As Double, Etmax As Double, Etmin As Double, Etmed As Double
    Dim Sen As Double, Senmax As Double, Senmin As Double, Senmed As Double
    ncmax = 10 ' Num. max. corpi
    lo = 639: lv = 479 ' Limiti grafici
    Dim um#(ncmax), ur#(ncmax), ud#(ncmax)
    Dim um$(ncmax), ur$(ncmax), ud$(ncmax)
    Dim M0(ncmax) As Double, M(ncmax) As Double
    Dim R0(ncmax) As Double, R(ncmax) As Double, Rg(ncmax) As Integer
    Dim x0(ncmax) As Double, y0(ncmax) As Double, z0(ncmax) As Double
    Dim x(ncmax) As Double, y(ncmax) As Double, z(ncmax) As Double
    Dim xa(ncmax) As Double, ya(ncmax) As Double, za(ncmax) As Double
    Dim xf(ncmax) As Double, yf(ncmax) As Double, zf(ncmax) As Double
    Dim xx1(ncmax) As Integer, yy1(ncmax) As Integer
    Dim xx2(ncmax) As Integer, yy2(ncmax) As Integer
    Dim zz1(ncmax) As Integer, zz2(ncmax) As Integer, prof(ncmax) As Integer
    Dim v(ncmax) As Double, vt(ncmax) As Double
    Dim vx(ncmax) As Double, vy(ncmax) As Double, vz(ncmax) As Double
    Dim dirxy(ncmax) As Double, dirz(ncmax) As Double, dz(ncmax) As Double
    Dim dmax0(ncmax) As Double, Dist(ncmax, ncmax) As Double, pdist(ncmax, ncmax) As Double
    Dim Apo(ncmax) As Double, Per(ncmax) As Double, SM(ncmax) As Double
    Dim AZ(ncmax) As Single, EL(ncmax) As Single
    Dim AN(ncmax) As Single
    Dim Ep(ncmax, ncmax) As Double, Ec(ncmax) As Double
    Dim Epg(lo) As Integer, Epmedg(lo) As Integer
    Dim Ecg(lo) As Integer, Ecmedg(lo) As Integer
    Dim Etg(lo) As Integer, Etmedg(lo) As Integer
    Dim Seng(lo) As Integer, Senmedg(lo) As Integer
    Randomize Timer
    
    '
    '     - Costanti -
    G = .0000000000667259# ' Costante gravitaz. universale
    VC = 299792458 ' Velocita' luce
    AL = VC * 31558149.5# ' Anno-luce
    pr = 100 ' Precisione
    pigreco = 3.14159265359#
    '
    ' parametri satellite geostazionario
    MT = 5.977D+24: RT = 6370000: Dsatg = 42170662
    ' parametri Sole-Terra
    MS = 1.989D+30: RS = 695980000: DTerra = 149598500000#
    '
    '     - Inizio -
    Screen 9, , 0: _FullScreen: L = 7: C = 1: nc = 2
    Cls: Locate 2, 20: Color 11: Print " NEWTON ": Color 2
    Print " - Simulazione del moto di molti corpi celesti"
    Print " " + String$(64, 196): Print Space$(46) + "(Premere un tasto)"
    Window (0, lv)-(lo, 0)
    GoSub comandi
    Cls: PCopy 0, 1: qq = 1
    '
    parametri:
    Cls: nc = 2: I1 = 0: I2 = 0
    x = 0: y = pigreco: z = 0 ' angoli di Eulero
    For N1 = 1 To nc: M0(N1) = 0: Next N1
    Color 5
    Locate 2, 65: Print "           ": Locate 2, 1
    Print "1) Satellite geostazionario "
    Print "2) Sole-Terra "
    Print "3) Random "
    Print "4) Inserisci parametri "
    100 sel$ = InKey$: If sel$ = "" Then 100
    If Asc(sel$) < 49 Or Asc(sel$) > 52 Then 100
    If sel$ = "1" Then
        nc = 2: M0(1) = MT: R0(1) = RT: M0(2) = 1: R0(2) = 1
        Mtot = M0(1) + M0(2)
        dmax0(2) = Dsatg
        Dist(1, 2) = dmax0(2): Dist(2, 1) = Dist(1, 2)
        GoSub velorb
        dirxy(1) = 0: dirxy(2) = 180: dirz(1) = 0: dirz(2) = 0
        dq0 = 1.5: Mmax = 1: Dmax = 2
        For N1 = 1 To nc: ud#(N1) = 1000: ud$(N1) = "Km": Next N1
        um#(1) = MT: ur#(1) = RT: um#(2) = 1: ur#(2) = 1
        um$(1) = "mt": ur$(1) = "rt": um$(2) = "Kg": ur$(2) = "metri"
        dq = dq0: GoSub pargraf ' par. grafici
        GoSub calcint ' calcola interv. tempo
        GoTo 110
    ElseIf sel$ = "2" Then
        nc = 2: M0(1) = MS: R0(1) = RS: M0(2) = MT: R0(2) = RT
        Mtot = M0(1) + M0(2)
        dmax0(2) = DTerra
        Dist(1, 2) = dmax0(2): Dist(2, 1) = Dist(1, 2)
        GoSub velorb
        dirxy(1) = 0: dirxy(2) = 180: dirz(1) = 0: dirz(2) = 0
        dq0 = 1.5: Mmax = 1: Dmax = 2
        For N1 = 1 To nc: ud#(N1) = 1000: ud$(N1) = "Km": Next N1
        um#(1) = MS: ur#(1) = RS: um#(2) = MT: ur#(2) = RT
        um$(1) = "ms": ur$(1) = "rs": um$(2) = "mt": ur$(2) = "rt"
        dq = dq0: GoSub pargraf: GoSub calcint
        GoTo 110
    ElseIf sel$ = "3" Then nc = 2: I1 = 1: GoTo 130
    Else GoTo 130
    End If
    110 Locate 2, 1: Print "massa  I astro    ="; M0(1) / um#(1); um$(1); "   "
    Print "raggio I astro    ="; R0(1) / ur#(1); ur$(1); "   "
    Print "massa  2 astro    ="; M0(2) / um#(2); um$(2)
    Print "raggio 2 astro    ="; R0(2) / ur#(2); ur$(2)
    Print "dim. quadro       ="; dq
    Print "distanza          ="; dmax0(2) / 1000; "Km   "
    Print "dir.xy I astro    = 180 gradi   "
    Print "dir.z  I astro    =   0 gradi   "
    Print "veloc. I astro    ="; v(1) / 1000; "Km/sec   "
    Print "dir.xy 2 astro    =   0 gradi   "
    Print "dir.z  2 astro    =   0 gradi   "
    Print "veloc. 2 astro    ="; v(2) / 1000; "Km/sec   "
    Print "interv.(sec)      ="; it; "sec   "
    Locate 16, 1: Print " Confermi (s) ? "
    120 A$ = InKey$: If A$ = "" Then 120
    If A$ <> "n" Then 190
    130 dq0 = 1
    140 Cls: PCopy 1, 0: Locate 2, 1
    L = 2: C = 27: Print "Numero astri ( max."; ncmax; "):"; nc
    GoSub scrivi: nc = N
    If nc < 2 Then nc = 2: GoTo 140
    If nc > ncmax Then nc = ncmax: GoTo 140
    If I1 Then dq = dq0: GoSub randomizza
    PCopy 1, 0: Locate 2, 1
    Print "Parametri (l'astro piu' massiccio diventa I astro)"
    Mtot = 0
    For N1 = 1 To nc
        GoSub unit1
        um#(N1) = um#: ur#(N1) = ur#: um$(N1) = um$: ur$(N1) = ur$
        M0(N1) = um#(N1): R0(N1) = ur#(N1)
        dmax0(N1) = 0: dirxy(N1) = 0: dirz(N1) = 0
        170 Locate 8, 1
        L = 8: C = 17: Print "massa"; N1; " ("; um$(N1); ") ="; M0(N1) / um#(N1);
        Print Space$(23 - Len(Str$(M0(N1) / um#(N1))))
        GoSub scrivi: M0(N1) = N * um#(N1)
        If M0(N1) <= 0 Then M0(N1) = um#(N1): GoTo 170
        M(N1) = M0(N1): Mtot = Mtot + M(N1)
        180 L = 9: C = 18: Print "raggio"; N1; " ("; ur$(N1); ") ="; R0(N1) / ur#(N1);
        Print Space$(23 - Len(Str$(R0(N1) / ur#(N1))))
        GoSub scrivi: R0(N1) = N * ur#(N1)
        If R0(N1) <= 0 Then R0(N1) = ur#(N1): Locate 7, 1: GoTo 180
        R(N1) = R0(N1)
        Lc = 11: GoSub angdir
    PCopy 1, 0: Next N1
    GoSub astromax: PCopy 1, 0
    If nc = 2 And M0(1) = MT And M0(2) < MT Then GoSub qudist Else GoSub distmax
    PCopy 1, 0
    For N1 = 1 To nc
        If nc = 2 Then GoSub velorb: GoTo 181
        v(N1) = Sqr(G * (M0(Mmax) ^ 2) / (M0(Mmax) + M0(N1)) / dmax0(Dmax))
        v(Mmax) = 0
        181 GoSub veloc
    Next N1
    dq = dq0: GoSub dimquad ' dim. quadro schermo
    GoSub interv ' precisa intervallo
    Locate 18, 1
    Print "Scegli assetto (n) ?"
    182 A$ = InKey$: If A$ = "" Then 182
    If A$ = "s" Then GoSub assetto: I1 = 1
    
    '
    '     - Inizializzazione -
    190 For N1 = 1 To nc
        xx1(N1) = -1: yy1(N1) = -1: xx2(N1) = -1: yy2(N1) = -1
    Next N1
    200 Screen 9, , 1: Cls
    210 Screen 9, , 0: PCopy 1, 0
    If I1 Then 240
    If I2 Then 260
    220 x0(Mmax) = 0: y0(Mmax) = 0: z0(Max) = 0
    For N1 = 1 To nc: If N1 = Mmax Then 230
        If N1 <> Mmax Then AZ(N1) = pigreco / 2: EL(N1) = 0
        x0(N1) = dmax0(N1) * Cos(AZ(N1)) * Cos(EL(N1))
        If x0(N1) < .1# Then x0(N1) = 0
        y0(N1) = dmax0(N1) * Sin(AZ(N1)) * Cos(EL(N1))
        If y0(N1) < .1# Then y0(N1) = 0
        z0(N1) = dmax0(N1) * Sin(EL(N1))
        If z0(N1) < .1# Then z0(N1) = 0
    230 Next N1
    240 If diag1 Or diag2 Or diag3 Then PCopy 1, 0
    If Nt = 1 Then Cls: PCopy 0, 1
    For N1 = 1 To nc
        M(N1) = M0(N1): R(N1) = R0(N1): Rg(N1) = R(N1) * kx
        dz(N1) = dirz(N1) * pigreco / 180
        vx(N1) = v(N1) * Cos(dz(N1)) * Cos(dirxy(N1) * pigreco / 180)
        vy(N1) = v(N1) * Cos(dz(N1)) * Sin(dirxy(N1) * pigreco / 180)
        vz(N1) = v(N1) * Sin(dz(N1))
        vt(N1) = Sqr(vx(N1) ^ 2 + vy(N1) ^ 2 + vz(N1) ^ 2)
        xx1(N1) = -1: yy1(N1) = -1: xx2(N1) = -1: yy2(N1) = -1
    Next N1: GoSub astromax
    For N1 = 1 To nc: If M(N1) = 0 Then 250
        x(N1) = x0(N1): y(N1) = y0(N1): z(N1) = z0(N1)
        Apo(N1) = 0: Per(N1) = dmax0(N1)
        If nc = 2 Then AN = AZ(N1)
        For N2 = 1 To nc: If N2 = N1 Or M(N2) = 0 Then 241
            Dist(N1, N2) = Sqr((x(N1) - x(N2)) ^ 2 + (y(N1) - y(N2)) ^ 2 + (z(N1) - z(N2)) ^ 2)
            Dist(N2, N1) = Dist(N1, N2)
            If Nt = 2 Then pdist(N1, N2) = Dist(N1, N2): pdist(N2, N1) = pdist(N1, N2)
        241 Next N2
    250 Next N1
    ao = 320.5: oo = 240.5 ' coord. centro schermo
    t = 0: I1 = 0: I2 = 0: Nt = nc
    pm = 0: primo = 0
    Epmax = 0: Epmin = -1D+300
    Ecmax = 0: Ecmin = 1D+300
    Etmax = 0: Etmin = -1D+300
    Senmax = 0: Senmin = -1D+300
    '
    '     - Calcolo matematico -
    260 vmax2 = 0
    For N1 = 1 To nc: If M0(N1) = 0 Then 261
        If vt(N1) > vmax2 Then vmax2 = vt(N1)
    261 Next N1
    ' - Calcolo precisione tempo -
    If vmax2 > vmax1 Then dt = Abs(it / pr * vmax1 / vmax2) Else dt = Abs(it / pr)
    For k = 1 To pr
        For N1 = 1 To nc: If M(N1) = 0 Then 280
            For N2 = 1 To nc: If N2 = N1 Or M(N2) = 0 Then 270
                Dist(N1, N2) = Sqr((x(N1) - x(N2)) ^ 2 + (y(N1) - y(N2)) ^ 2 + (z(N1) - z(N2)) ^ 2)
                Dist(N2, N1) = Dist(N1, N2)
                If Nt = 2 Then pdist(N1, N2) = Dist(N1, N2): pdist(N2, N1) = pdist(N1, N2)
            270 Next N2
        280 Next N1
        For N1 = 1 To nc: If M(N1) = 0 Then 291
            xa(N1) = 0: ya(N1) = 0: za(N1) = 0
            If Nt > 1 Then
                For N2 = 1 To nc: If N2 = N1 Or M(N2) = 0 Then 290
                    acc = G * M(N2) / Dist(N1, N2) ^ 3
                    xa(N1) = xa(N1) + acc * (x(N2) - x(N1))
                    ya(N1) = ya(N1) + acc * (y(N2) - y(N1))
                    za(N1) = za(N1) + acc * (z(N2) - z(N1))
                290 Next N2
            End If
        291 Next N1
        For N1 = 1 To nc: If M(N1) = 0 Then 292
            If Nt > 1 Then
                vx(N1) = vx(N1) + xa(N1) * dt
                vy(N1) = vy(N1) + ya(N1) * dt
                vz(N1) = vz(N1) + za(N1) * dt
                vt(N1) = Sqr(vx(N1) ^ 2 + vy(N1) ^ 2 + vz(N1) ^ 2)
            End If
        292 Next N1
        For N1 = 1 To nc: If M(N1) = 0 Then 300
            xf(N1) = x(N1) + vx(N1) * dt
            yf(N1) = y(N1) + vy(N1) * dt
            zf(N1) = z(N1) + vz(N1) * dt
        300 Next N1
        For N1 = 1 To nc: If M(N1) = 0 Then 330
            x(N1) = xf(N1): y(N1) = yf(N1): z(N1) = zf(N1)
            If nc = 2 Then
                AP = AN: If x(2) = 0 Then AN = pigreco / 2 Else AN = Atn(y(2) / x(2))
            End If
            '
            '     - Controlli -
            For N2 = 1 To nc: If N2 = N1 Or M(N2) = 0 Then 320
                Dist(N1, N2) = Sqr((x(N1) - x(N2)) ^ 2 + (y(N1) - y(N2)) ^ 2 + (z(N1) - z(N2)) ^ 2)
                Dist(N2, N1) = Dist(N1, N2)
                If N1 = Mmax Then 310
                fndistmax = Sqr((x(N1) - x(Mmax)) ^ 2 + (y(N1) - y(Mmax)) ^ 2 + (z(N1) - z(Mmax)) ^ 2)
                If fndistmax > Apo(N1) Then Apo(N1) = fndistmax
                If fndistmax < Per(N1) Then Per(N1) = fndistmax
                SM(N1) = (Apo(N1) + Per(N1)) / 2
                310 If Dist(N1, N2) > R(N1) + R(N2) Then 320
                '
                ' altrimenti...
                '
                '     - Collisione-fusione -
                Locate 2, 65: Color 15
                Sound 200, .1: Print "C R A S H !"
                Locate 4, 1: Color 3
                Print "Tempo =";
                tanni = 0: tmesi = 0: tgiorni = 0: tore = 0: tminuti = 0: tsecondi = 0
                tanni = Int(t / 31536000)
                tmesi = Int((t - (tanni * 31536000)) / 2592000)
                tgiorni = Int((t - (tanni * 31536000) - (tmesi * 2592000)) / 86400)
                tore = Int((t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400)) / 3600)
                tminuti = Int((t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400) - (tore * 3600)) / 60)
                tsecondi = Int(t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400) - (tore * 3600) - (tminuti * 60))
                Print tanni; "yr"; tmesi; "mo"; tgiorni; "d"; tore; "h"; tminuti; "m"; tsecondi; "s ";
                Print "="; t; "sec "
                Print "Num. corpi ="; Nt
                Print: Print "Massa 1 ="; M(N1); "Kg"
                Print "Vt 1 =";
                If vt(N1) > 1000 Then Print vt(N1) / 1000; "Km/s " Else Print vt(N1); "m/s "
                Print: Print "Massa 2 ="; M(N2); "Kg"
                Print "Vt 2 =";
                If vt(N2) > 1000 Then Print vt(N2) / 1000; "Km/s " Else Print vt(N2); "m/s "
                If Nt = 2 Then
                    Print: Print "En. cin. scontro = ";
                    Ecs = (M(N1) * vt(N1) ^ 2 + M(N2) * vt(N2) ^ 2) / 2: Print Ecs; "Joule ";
                End If
                Do: Loop While InKey$ = ""
                Cls: PCopy 0, 1: Nt = Nt - 1
                m1 = N1: m2 = N2: If M(N2) > M(N1) Then Swap m1, m2
                den = M(m1) * 3 / 4 / pigreco / R(m1) ^ 3
                vx(m1) = (M(m1) * vx(m1) + M(m2) * vx(m2)) / (M(m1) + M(m2)): vx(m2) = 0
                vy(m1) = (M(m1) * vy(m1) + M(m2) * vy(m2)) / (M(m1) + M(m2)): vy(m2) = 0
                vz(m1) = (M(m1) * vz(m1) + M(m2) * vz(m2)) / (M(m1) + M(m2)): vz(m2) = 0
                M(m1) = M(m1) + M(m2): M(m2) = 0
                R(m1) = ((M(m1) / den) * 3 / 4 / pigreco) ^ (1 / 3)
                Rg(m1) = R(m1) * kx
                GoSub astromax: pm = 0: primo = 0
                Epmax = 0: Epmin = -1D+300
                Ecmax = 0: Ecmin = 1D+300
                Etmax = 0: Etmin = -1D+300
                Senmax = 0: Senmin = -1D+300
                GoTo 260
                '
            320 Next N2
        330 Next N1
    Next k: t = t + it * vmax2 / vmax1
    '
    ' calcolo energia istantanea, massima, minima e media
    Ep = 0: Ec = 0
    For N1 = 1 To nc: If M(N1) = 0 Then 350
        For N2 = 1 To nc: If N2 = N1 Or M(N2) = 0 Then 340
            fndist = Sqr((x(N1) - x(N2)) ^ 2 + (y(N1) - y(N2)) ^ 2 + (z(N1) - z(N2)) ^ 2)
            Ep(N1, N2) = -G * M(N1) * M(N2) / fndist
            Ep = Ep + Ep(N1, N2)
        340 Next N2
        Ec(N1) = vt(N1) ^ 2 * M(N1) / 2: Ec = Ec + Ec(N1)
    350 Next N1
    Ep = Ep / 2: Etot = Ep + Ec
    If Abs(Ep) > Abs(Epmax) Then Epmax = Ep
    If Abs(Ep) < Abs(Epmin) Then Epmin = Ep
    If Ec > Ecmax Then Ecmax = Ec
    If Ec < Ecmin Then Ecmin = Ec
    If Abs(Etot) > Abs(Etmax) Then Etmax = Etot
    If Abs(Etot) < Abs(Etmin) Then Etmin = Etot
    Epmed = (Epmax + Epmin) / 2
    Ecmed = (Ecmax + Ecmin) / 2
    Etmed = (Etmax + Etmin) / 2
    Sen = Etot / Ec ' entropia
    If Abs(Sen) > Abs(Senmax) Then Senmax = Sen
    If Abs(Sen) < Abs(Senmin) Then Senmin = Sen
    Senmed = (Senmax + Senmin) / 2
    '
    ' fsd = fattore scala diagrammi
    If primo = 0 Then
        If Epmin > -1D+300 Then
            fsd1 = (lv - 1) / Abs(Epmin) / 5
        Else fsd1 = (lv - 1) / Abs(Ep) / 5
        End If
        If Ep = 0 Then fsd1 = (lv - 1) / Ec / 2
        If Senmin > -1D-300 Then
            fsd2 = 0.9 * (lv - 1) / Abs(Senmin) / 2
        Else fsd2 = 0.9 * (lv - 1) / Abs(Sen) / 2
        End If
        If Sen = 0 Then fsd2 = 0.9 * (lv - 1) / 2
        primo = 1
    End If
    Epg = Ep * fsd1 + (lv - 1) / 2
    If Epg > lv - 1 Or Epg < 0 Then fsd1 = fsd1 / 2: pm = lo
    Epmedg = Epmed * fsd1 + (lv - 1) / 2
    If Epmedg > lv - 1 Or Epmedg < 0 Then fsd1 = fsd1 / 2: pm = lo
    Ecg = Ec * fsd1 + (lv - 1) / 2
    If Ecg > lv - 1 Or Ecg < 0 Then fsd1 = fsd1 / 2: pm = lo
    Ecmedg = Ecmed * fsd1 + (lv - 1) / 2
    If Ecmedg > lv - 1 Or Ecmedg < 0 Then fsd1 = fsd1 / 2: pm = lo
    Etg = Etot * fsd1 + (lv - 1) / 2
    If Etg > lv - 1 Or Etg < 0 Then fsd1 = fsd1 / 2: pm = lo
    Etmedg = Etmed * fsd1 + (lv - 1) / 2
    If Etmedg > lv - 1 Or Etmedg < 0 Then fsd1 = fsd1 / 2: pm = lo
    Seng = Sen * fsd2 + (lv - 1) / 2
    If Seng > lv - 1 Or Seng < 0 Then fsd2 = fsd2 / 2: pm = lo
    Senmedg = Senmed * fsd2 + (lv - 1) / 2
    If Senmedg > lv - 1 Or Senmedg < 0 Then fsd2 = fsd2 / 2: pm = lo
    pm = pm + 1
    If pm > lo / 6 Then
        pm = 0: If diag1 Or diag2 Or diag3 Then Cls: PCopy 1, 0
    End If
    Epg(pm) = Epg: Epmedg(pm) = Epmedg
    Ecg(pm) = Ecg: Ecmedg(pm) = Ecmedg
    Etg(pm) = Etg: Etmedg(pm) = Etmedg
    Seng(pm) = Seng: Senmedg(pm) = Senmedg
    '
    GoSub grafica ' <- routine grafica
    '
    A$ = InKey$: If A$ = "" Then 260
    If A$ = "e" Then L = 3: C = 3: GoSub comandi: I2 = 1: GoTo 210
    If A$ = "p" Then GoTo parametri
    If A$ = "q" Then GoSub dimquad: GoTo 190
    If A$ = "d" Then
        If nc = 2 And M0(2) < MT Then GoSub qudist Else GoSub distmax
        dq = dq0: GoSub pargraf: GoTo 190
    End If
    If A$ = "v" Then
        For N1 = 1 To nc: GoSub veloc: Next N1
        dq = dq0: GoSub pargraf: I1 = 1: GoTo 190
    End If
    If A$ = "g" Then
        I2 = 1: For N1 = 1 To nc: Lc = 21: GoSub angdir: Next N1
        dq = dq0: GoSub pargraf: I1 = 1: GoTo 190
    End If
    If A$ = "a" Then I2 = 1: GoSub assetto: I1 = 1: GoTo 190
    If A$ = "i" Then GoSub interv: I2 = 1: GoTo 190
    If A$ = "t" Then
        For N1 = 1 To nc: vx(N1) = -vx(N1): vy(N1) = -vy(N1): Next N1
        it = -it
    End If
    If A$ = "r" Then I1 = 1: GoTo 130
    If A$ = "n" Then dq = dq0: GoSub pargraf: I1 = 1: GoTo 200
    If A$ = "c" Then Cls: PCopy 0, 1
    If A$ = "h" Then GoSub 4000: GoTo 260
    If A$ = "<" Then
        dq = dq * 2
        GoSub pargraf: I2 = 1: GoTo 190
    End If
    If A$ = ">" Then
        dq = dq / 2: If dq < .001 Then dq = .001
        GoSub pargraf: I2 = 1: GoTo 190
    End If
    If Asc(A$) > 47 And Asc(A$) < 58 Then GoSub dati: I2 = 1: GoTo 210
    If A$ = " " Then GoSub globali: I2 = 1: GoTo 210
    If A$ = "!" Then
        diag1 = -(diag1 = 0): diag2 = 0
        pm = lo
        Cls: PCopy 1, 0
    End If
    If Asc(A$) = 34 Then
        diag2 = -(diag2 = 0): diag1 = 0
        pm = lo
        Cls: PCopy 1, 0
    End If
    If A$ = "�" Then
        diag3 = -(diag3 = 0): diag1 = 0: diag2 = 0
        pm = lo
        Cls: PCopy 1, 0
    End If
    If A$ = Chr$(0) + Chr$(59) Or A$ = Chr$(0) + Chr$(60) Or A$ = Chr$(0) + Chr$(61) Or A$ = Chr$(0) + Chr$(62) Then
        If A$ = Chr$(0) + Chr$(59) Then qq = 1
        If A$ = Chr$(0) + Chr$(60) Then qq = 2
        If A$ = Chr$(0) + Chr$(61) Then qq = 3
        dq = dq0: GoSub pargraf: I2 = 1: GoTo 190
    End If
    If A$ = "x" Or A$ = "X" Or A$ = "y" Or A$ = "Y" Or A$ = "z" Or A$ = "Z" Then
        If A$ = "x" Then x = x + 5 * pigreco / 180
        If x > (355 * pigreco / 180) Then x = 0
        If A$ = "X" Then x = x - 5 * pigreco / 180
        If x < -(355 * pigreco / 180) Then x = 0
        If A$ = "y" Then y = y + 5 * pigreco / 180
        If y > (445 * pigreco / 180) Then y = pigreco / 2
        If A$ = "Y" Then y = y - 5 * pigreco / 180
        If y < -(265 * pigreco / 180) Then y = pigreco / 2
        If A$ = "z" Then z = z - 5 * pigreco / 180
        If z < -(355 * pigreco / 180) Then z = 0
        If A$ = "Z" Then z = z + 5 * pigreco / 180
        If z > (355 * pigreco / 180) Then z = 0
        dq = dq0: GoSub pargraf: I2 = 1: GoTo 190
    End If
    If Len(A$) >= 2 Then
        AK = Asc(Right$(A$, 1))
        If AK = 72 Then oo = oo + 10
        If AK = 80 Then oo = oo - 10
        If AK = 75 Then ao = ao - 10
        If AK = 77 Then ao = ao + 10
        dq = dq0: GoSub pargraf: I2 = 1: GoTo 190
    End If
    If A$ = Chr$(27) Then Screen 0: End
    GoTo 260
    '
    '     - Routines -
    grafica:
    
    ' vista panoramica
    If qq = 1 Then
        Screen 9, , 0
        Locate 2, 1: Color 15: Print "Vista panoramica"
        For N1 = 1 To nc: If M(N1) = 0 Then 1001
            If xx1(N1) > 0 And xx1(N1) < lo And yy1(N1) > 0 And yy1(N1) < lv Then
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    Circle (xx1(N1), yy1(N1)), Rg(N1), 0, , , 1
                    If nc = 2 Then
                        If Sqr((xx1(1) - xx1(2)) ^ 2 + (yy1(1) - yy1(2)) ^ 2) > Rg(1) + Rg(2) Then
                            PSet (xx1(N1), yy1(N1)), 5
                            GoTo 1000
                        End If
                        If zz1(1) < zz1(2) Then
                            PSet (xx1(1), yy1(1)), 5
                        Else PSet (xx1(2), yy1(2)), 5
                        End If
                    End If
                1000 Next scr
            End If
        1001 Next N1
        1002 For N1 = 1 To nc: If M(N1) = 0 Then 1003
            xx1(Mmax) = ((Cos(x) * Cos(z) + Sin(x) * Cos(y) * Sin(z)) * (x(Mmax) * 0.75) + (-Sin(x) * Cos(z) + Cos(x) * Cos(y) * Sin(z)) * y(Mmax) + (Sin(y) * Sin(z)) * z(Mmax)) * kx + ao
            yy1(Mmax) = ((-Cos(x) * Sin(z) - Sin(x) * Cos(y) * Cos(z)) * (x(Mmax) * 0.75) + (Sin(x) * Sin(z) - Cos(x) * Cos(y) * Cos(z)) * y(Mmax) + (Sin(y) * Cos(z)) * z(Mmax)) * kx + oo
            xx1(N1) = ((Cos(x) * Cos(z) + Sin(x) * Cos(y) * Sin(z)) * (x(N1) * 0.75) + (-Sin(x) * Cos(z) + Cos(x) * Cos(y) * Sin(z)) * y(N1) + (Sin(y) * Sin(z)) * z(N1)) * kx + ao
            yy1(N1) = ((-Cos(x) * Sin(z) - Sin(x) * Cos(y) * Cos(z)) * (x(N1) * 0.75) + (Sin(x) * Sin(z) - Cos(x) * Cos(y) * Cos(z)) * y(N1) + (Sin(y) * Cos(z)) * z(N1)) * kx + oo
            zz1(N1) = ((Sin(x) * Sin(y)) * (x(N1) * 0.75) + (Cos(x) * Sin(y)) * y(N1) + Cos(y) * z(N1)) * kx + ao
            If xx1(Mmax) < -1 Or xx1(Mmax) > lo Or yy1(Mmax) < -1 Or yy1(Mmax) > lv Then
                If xx1(Mmax) < -1 Then ao = ao + lo / 2
                If xx1(Mmax) > lo Then ao = ao - lo / 2
                If yy1(Mmax) < -1 Then oo = oo + lv / 2
                If yy1(Mmax) > lv Then oo = oo - lv / 2
                Screen 9, , 0, 0: Cls: PCopy 0, 1: GoTo 1002
            End If
        1003 Next N1
        GoSub zprof1
        For N1 = 1 To nc: If M(prof(N1)) = 0 Then 1006
            If xx1(N1) < -1 Or xx1(N1) > lo Or yy1(N1) < -1 Or yy1(N1) > lv Then
                dq = dq * 2
                GoSub pargraf: I2 = 1: GoTo 190
            End If
            For N2 = 2 To nc: If M(prof(N2)) = 0 Then 1005
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    col = 15: If nc > 2 And nc < 14 Then col = N1 + 1
                    If M(prof(N1)) >= MT Then col = 9
                    If M(prof(N1)) >= MS Then col = 14
                    If M(prof(N2)) >= MT Then col = 9
                    If M(prof(N2)) >= MS Then col = 14
                    If N1 = Mmax And sel$ = "2" Then col = 14
                    If N1 = Mmax And sel$ = "3" Then col = 15
                    If Sqr((xx1(prof(N1)) - xx1(prof(N2))) ^ 2 + (yy1(prof(N1)) - yy1(prof(N2))) ^ 2) > Rg(prof(N1)) + Rg(prof(N2)) Then
                        Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                        GoTo 1004
                    End If
                    If zz1(prof(N1)) < zz1(prof(N2)) Then
                        If Rg(prof(N2)) >= Rg(prof(N1)) Then
                            Circle (xx1(prof(N2)), yy1(prof(N2))), Rg(prof(N2)), col, , , 1
                        End If
                        Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                    Else
                        If Rg(prof(N1)) >= Rg(prof(N2)) Then
                            Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                        End If
                        Circle (xx1(prof(N2)), yy1(prof(N2))), Rg(prof(N2)), col, , , 1
                    End If
                1004 Next scr
            1005 Next N2
        1006 Next N1
        Screen 9, , 0
        Locate 21, 69: Color 15: Print "x ="; CInt(x * 180 / pigreco);
        Locate 22, 69: Color 15: Print "y ="; CInt(y * 180 / pigreco - 90);
        Locate 23, 69: Color 15: Print "z ="; CInt(z * 180 / pigreco);
    End If
    
    
    ' assonometria
    If qq = 2 Then
        Screen 9, , 0
        Locate 2, 1: Color 15: Print "Assonometria"
        Locate 24: Print "Y";
        Locate 19, 80: Print "X"
        Locate 1, 18: Print "Z"
        Line (0, 0)-(150, 150), 8
        Line (150, 150)-(150, lv), 8
        Line (150, 150)-(lo, 150), 8
        For N1 = 1 To nc: If M(N1) = 0 Then 1041
            If xx1(N1) > 0 And xx1(N1) < lo And yy1(N1) > 0 And yy1(N1) < lv Then
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    Circle (xx1(N1), yy1(N1)), Rg(N1), 0, , , 1
                    If nc = 2 Then
                        If Sqr((xx1(1) - xx1(2)) ^ 2 + (yy1(1) - yy1(2)) ^ 2) > Rg(1) + Rg(2) Then
                            PSet (xx1(N1), yy1(N1)), 5
                            GoTo 1040
                        End If
                        If zz1(1) < zz1(2) Then
                            PSet (xx1(1), yy1(1)), 5
                        Else PSet (xx1(2), yy1(2)), 5
                        End If
                    End If
                1040 Next scr
            End If
        1041 Next N1
        1042 For N1 = 1 To nc: If M(N1) = 0 Then 1043
            xx1(Mmax) = ((x(Mmax) * 0.75) - (y(Mmax) * .25)) * kx + ao
            yy1(Mmax) = (z(Mmax) + (y(Mmax) * .25)) * kx + oo
            xx1(N1) = ((x(N1) * 0.75) - (y(N1) * .25)) * kx + ao
            yy1(N1) = (z(N1) + (y(N1) * .25)) * kx + oo
            zz1(N1) = y(N1) * kx + ao
            If xx1(Mmax) < -1 Or xx1(Mmax) > lo Or yy1(Mmax) < -1 Or yy1(Mmax) > lv Then
                If xx1(Mmax) < -1 Then ao = ao + lo / 2
                If xx1(Mmax) > lo Then ao = ao - lo / 2
                If yy1(Mmax) < -1 Then oo = oo + lv / 2
                If yy1(Mmax) > lv Then oo = oo - lv / 2
                Screen 9, , 0, 0: Cls: PCopy 0, 1: GoTo 1042
            End If
        1043 Next N1
        GoSub zprof1
        For N1 = 1 To nc: If M(prof(N1)) = 0 Then 1046
            If xx1(N1) < -1 Or xx1(N1) > lo Or yy1(N1) < -1 Or yy1(N1) > lv Then
                dq = dq * 2
                GoSub pargraf: I2 = 1: GoTo 190
            End If
            For N2 = 2 To nc: If M(prof(N2)) = 0 Then 1045
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    col = 15: If nc > 2 And nc < 14 Then col = N1 + 1
                    If M(prof(N1)) >= MT Then col = 9
                    If M(prof(N1)) >= MS Then col = 14
                    If M(prof(N2)) >= MT Then col = 9
                    If M(prof(N2)) >= MS Then col = 14
                    If N1 = Mmax And sel$ = "2" Then col = 14
                    If N1 = Mmax And sel$ = "3" Then col = 15
                    If Sqr((xx1(prof(N1)) - xx1(prof(N2))) ^ 2 + (yy1(prof(N1)) - yy1(prof(N2))) ^ 2) > Rg(prof(N1)) + Rg(prof(N2)) Then
                        Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                        GoTo 1044
                    End If
                    If zz1(prof(N1)) < zz1(prof(N2)) Then
                        If Rg(prof(N2)) >= Rg(prof(N1)) Then
                            Circle (xx1(prof(N2)), yy1(prof(N2))), Rg(prof(N2)), col, , , 1
                        End If
                        Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                    Else
                        If Rg(prof(N1)) >= Rg(prof(N2)) Then
                            Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                        End If
                        Circle (xx1(prof(N2)), yy1(prof(N2))), Rg(prof(N2)), col, , , 1
                    End If
                1044 Next scr
            1045 Next N2
        1046 Next N1
    End If
    
    
    ' proiezioni x-y e x-z
    If qq = 3 Then
        Screen 9, , 0
        Locate 2, 2: Color 15: Print "Proiezione x-y e x-z"
        Color 5
        Locate 1, 2: Print "Y"
        Locate 25, 36: Print "X";
        Locate 1, 42: Print "Z"
        Locate 25, 80: Print "X";
        Line (1, lv)-(1, 1), 8
        Line -(280, 1), 8
        Line (320, lv)-(320, 1), 8
        Line -(630, 1), 8
        For N1 = 1 To nc: If M(N1) = 0 Then 1011
            If xx1(N1) > 0 And xx1(N1) < 315 And yy1(N1) > 0 And yy1(N1) < lv Then
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    Circle (xx1(N1), yy1(N1)), Rg(N1), 0, , , 1
                    If nc = 2 Then
                        If Sqr((xx1(1) - xx1(2)) ^ 2 + (yy1(1) - yy1(2)) ^ 2) > Rg(1) + Rg(2) Then
                            PSet (xx1(N1), yy1(N1)), 5
                            GoTo 1010
                        End If
                        If zz1(1) < zz1(2) Then
                            PSet (xx1(1), yy1(1)), 5
                        Else PSet (xx1(2), yy1(2)), 5
                        End If
                    End If
                1010 Next scr
            End If
        1011 Next N1
        For N1 = 1 To nc: If M(N1) = 0 Then 1013
            If xx2(N1) > 320 And xx2(N1) < lo And yy2(N1) > 0 And yy2(N1) < lv Then
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    Circle (xx2(N1), yy2(N1)), Rg(N1), 0, , , 1
                    If nc = 2 Then
                        If Sqr((xx2(1) - xx2(2)) ^ 2 + (yy2(1) - yy2(2)) ^ 2) > Rg(1) + Rg(2) Then
                            PSet (xx2(N1), yy2(N1)), 5
                            GoTo 1012
                        End If
                        If zz2(1) < zz2(2) Then
                            PSet (xx2(1), yy2(1)), 5
                        Else PSet (xx2(2), yy2(2)), 5
                        End If
                    End If
                1012 Next scr
            End If
        1013 Next N1
        1014 For N1 = 1 To nc: If M(N1) = 0 Then 1015
            xx1(Mmax) = (x(Mmax) * 0.75) * kx + ao / 2
            yy1(Mmax) = y(Mmax) * kx + oo
            xx1(N1) = (x(N1) * 0.75) * kx + ao / 2
            yy1(N1) = y(N1) * kx + oo
            zz1(N1) = z(N1) * kx + ao
            If xx1(Mmax) < -1 Or xx1(Mmax) > 315 Or yy1(Mmax) < -1 Or yy1(Mmax) > lv Then
                If xx1(Mmax) < -1 Then ao = ao + lo / 4
                If xx1(Mmax) > 315 Then ao = ao - lo / 4
                If yy1(Mmax) < -1 Then oo = oo + lv / 2
                If yy1(Mmax) > lv Then oo = oo - lv / 2
                Screen 9, , 0, 0: Cls: PCopy 0, 1: GoTo 1015
            End If
        1015 Next N1
        GoSub zprof1
        For N1 = 1 To nc: If M(prof(N1)) = 0 Then 1018
            If xx1(N1) < -1 Or xx1(N1) > 315 Or yy1(N1) < -1 Or yy1(N1) > lv Then
                dq = dq * 2
                GoSub pargraf: I2 = 1: GoTo 190
            End If
            For N2 = 2 To nc: If M(prof(N2)) = 0 Then 1017
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    col = 15: If nc > 2 And nc < 14 Then col = N1 + 1
                    If M(prof(N1)) >= MT Then col = 9
                    If M(prof(N1)) >= MS Then col = 14
                    If M(prof(N2)) >= MT Then col = 9
                    If M(prof(N2)) >= MS Then col = 14
                    If N1 = Mmax And sel$ = "2" Then col = 14
                    If N1 = Mmax And sel$ = "3" Then col = 15
                    If Sqr((xx1(prof(N1)) - xx1(prof(N2))) ^ 2 + (yy1(prof(N1)) - yy1(prof(N2))) ^ 2) > Rg(prof(N1)) + Rg(prof(N2)) Then
                        Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                        GoTo 1016
                    End If
                    If zz1(prof(N1)) < zz1(prof(N2)) Then
                        If Rg(prof(N2)) >= Rg(prof(N1)) Then
                            Circle (xx1(prof(N2)), yy1(prof(N2))), Rg(prof(N2)), col, , , 1
                        End If
                        Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                    Else
                        If Rg(prof(N1)) >= Rg(prof(N2)) Then
                            Circle (xx1(prof(N1)), yy1(prof(N1))), Rg(prof(N1)), col, , , 1
                        End If
                        Circle (xx1(prof(N2)), yy1(prof(N2))), Rg(prof(N2)), col, , , 1
                    End If
                1016 Next scr
            1017 Next N2
        1018 Next N1
        1019 For N1 = 1 To nc: If M(N1) = 0 Then 1020
            xx2(Mmax) = (x(Mmax) * 0.75) * kx + ao * 1.5
            yy2(Mmax) = z(Mmax) * kx + oo
            xx2(N1) = (x(N1) * 0.75) * kx + ao * 1.5
            yy2(N1) = z(N1) * kx + oo
            zz2(N1) = y(N1) * kx + ao
            If xx2(Mmax) < 320 Or xx2(Mmax) > lo Or yy2(Mmax) < -1 Or yy2(Mmax) > lv Then
                If xx2(Mmax) < 320 Then ao = ao + lo / 4
                If xx2(Mmax) > lo Then ao = ao - lo / 4
                If yy2(Mmax) < -1 Then oo = oo + lv / 2
                If yy2(Mmax) > lv Then oo = oo - lv / 2
                Screen 9, , 0, 0: Cls: PCopy 0, 1: GoTo 1020
            End If
        1020 Next N1
        GoSub zprof2
        For N1 = 1 To nc: If M(prof(N1)) = 0 Then 1023
            If xx2(prof(N1)) < 320 Or xx2(prof(N1)) > lo Or yy2(prof(N1)) < -1 Or yy2(prof(N1)) > lv Then
                dq = dq * 2
                GoSub pargraf: I2 = 1: GoTo 190
            End If
            For N2 = 2 To nc: If M(prof(N2)) = 0 Then 1022
                For scr = 0 To 1
                    Screen 9, , scr, 0
                    col = 15: If nc > 2 And nc < 14 Then col = N1 + 1
                    If M(prof(N1)) >= MT Then col = 9
                    If M(prof(N1)) >= MS Then col = 14
                    If M(prof(N2)) >= MT Then col = 9
                    If M(prof(N2)) >= MS Then col = 14
                    If N1 = Mmax And sel$ = "2" Then col = 14
                    If N1 = Mmax And sel$ = "3" Then col = 15
                    If Sqr((xx2(prof(N1)) - xx2(prof(N2))) ^ 2 + (yy2(prof(N1)) - yy2(prof(N2))) ^ 2) > Rg(prof(N1)) + Rg(prof(N2)) Then
                        Circle (xx2(prof(N1)), yy2(prof(N1))), Rg(prof(N1)), col, , , 1
                        GoTo 1021
                    End If
                    If zz2(prof(N1)) < zz2(prof(N2)) Then
                        If Rg(prof(N2)) >= Rg(prof(N1)) Then
                            Circle (xx2(prof(N2)), yy2(prof(N2))), Rg(prof(N2)), col, , , 1
                        End If
                        Circle (xx2(prof(N1)), yy2(prof(N1))), Rg(prof(N1)), col, , , 1
                    Else
                        If Rg(prof(N1)) >= Rg(prof(N2)) Then
                            Circle (xx2(prof(N1)), yy2(prof(N1))), Rg(prof(N1)), col, , , 1
                        End If
                        Circle (xx2(prof(N2)), yy2(prof(N2))), Rg(prof(N2)), col, , , 1
                    End If
                1021 Next scr
            1022 Next N2
        1023 Next N1
    End If
    
    
    Screen 9, , 0: Locate 24, 69: Color 15: Print "dq ="; dq;
    
    
    '             traccia ENERGIA TOTALE
    If diag1 Then
        Locate 25, 2: Color 15: Print "energia ";
        Color 14: Print "potenziale";: Color 15: Print ", ";
        Color 4: Print "cinetica";: Color 15: Print ", ";
        Color 2: Print "totale";
        Line (0, (lv - 1) / 2)-(lo, (lv - 1) / 2), 15
        ' en. potenziale istantanea
        PSet (0, Epg(0)), 14: fl% = 1
        For k = 1 To pm
            xx = k * 6
            If fl% Then PSet (xx, Epg(k)), 14: fl% = 0
            Line -(xx, Epg(k)), 14
        Next
        ' en. cinetica istantanea
        PSet (0, Ecg(0)), 4: fl% = 1
        For k = 1 To pm
            xx = k * 6
            If fl% Then PSet (xx, Ecg(k)), 4: fl% = 0
            Line -(xx, Ecg(k)), 4
        Next
        ' en. totale istantanea
        PSet (0, Etg(0)), 2: fl% = 1
        For k = 1 To pm
            xx = k * 6
            If fl% Then PSet (xx, Etg(k)), 2: fl% = 0
            Line -(xx, Etg(k)), 2
        Next
        ' en. potenziale media
        PSet (0, Epmedg(0)), 15
        For k = 1 To pm
            xx = k * 6
            PSet (xx, Epmedg(k)), 15
        Next
        ' en. cinetica media
        PSet (0, Ecmedg(0)), 15
        For k = 1 To pm
            xx = k * 6
            PSet (xx, Ecmedg(k)), 15
        Next
        ' en. totale media
        PSet (0, Etmedg(0)), 15
        For k = 1 To pm
            xx = k * 6
            PSet (xx, Etmedg(k)), 15
        Next
    End If
    '
    '             traccia ENTROPIA
    If diag2 Then
        Locate 25, 2: Color 9: Print "entropia (Etot/Ecin)";
        Line (0, (lv - 1) / 2)-(lo, (lv - 1) / 2), 15
        ' entropia istantanea
        PSet (0, Seng(0)), 9: fl% = 1
        For k = 1 To pm
            xx = k * 6
            If fl% Then PSet (xx, Seng(k)), 9: fl% = 0
            Line -(xx, Seng(k)), 9
        Next
        ' entropia media
        PSet (0, Senmedg(0)), 15
        For k = 1 To pm
            xx = k * 6
            PSet (xx, Senmedg(k)), 15
        Next
    End If
    '
    Return
    '
    dati:
    '     - Differenza angolare -
    If nc = 2 Then
        DA = Abs(AN - AP): If AN * AP > 0 Then 2000
        If AP < 0 Then AP = AP + pigreco
        DA = AP - AN
    End If
    '     - Stampa dati -
    2000 N1 = Val(A$): If N1 = 0 Then N1 = 10
    Locate 3, 1
    col = 5: If nc > 2 And nc < 14 Then col = N1 + 1
    If M(N1) >= MT Then col = 9
    If M(N1) >= MS Then col = 14
    If N1 = Mmax And sel$ = "3" Then col = 15
    Color col
    If M(N1) = 0 Then 2040
    Print "Tempo =";
    tanni = 0: tmesi = 0: tgiorni = 0: tore = 0: tminuti = 0: tsecondi = 0
    tanni = Int(t / 31536000)
    tmesi = Int((t - (tanni * 31536000)) / 2592000)
    tgiorni = Int((t - (tanni * 31536000) - (tmesi * 2592000)) / 86400)
    tore = Int((t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400)) / 3600)
    tminuti = Int((t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400) - (tore * 3600)) / 60)
    tsecondi = Int(t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400) - (tore * 3600) - (tminuti * 60))
    Print tanni; "yr"; tmesi; "mo"; tgiorni; "d"; tore; "h"; tminuti; "m"; tsecondi; "s ";
    Print "="; t; "sec "
    Print "Num. corpi ="; Nt
    Print N1; ") Massa ="; M(N1); "Kg"
    If N1 = Mmax Or M(Mmax) = 0 Then 2020
    Print "     Dist.I astro =";: fndistmax = Sqr((x(N1) - x(Mmax)) ^ 2 + (y(N1) - y(Mmax)) ^ 2 + (z(N1) - z(Mmax)) ^ 2)
    If fndistmax < 1000 Then Print fndistmax; "m"
    If fndistmax < AL Then Print fndistmax / 1000; "Km" Else Print fndistmax / AL; "al"
    Print "     Apoastro     =";
    If Apo(N1) < 1000 Then Print Apo(N1); "m"
    If Apo(N1) < AL Then Print Apo(N1) / 1000; "Km" Else Print Apo(N1) / AL; "al"
    Print "     Periastro    =";
    If Per(N1) < 1000 Then Print Per(N1); "m"
    If Per(N1) < AL Then Print Per(N1) / 1000; "Km" Else Print Per(N1) / AL; "al"
    If Nt > 2 Or M(Mmax) <> MT Then 2010
    Print "     Quota        =";: Quota = fndistmax - R(Mmax)
    If Quota >= 1000 Then Print Quota / 1000; "Km" Else Print Quota; "m"
    Print "     QApo         =";: QApo = Apo(N1) - R(Mmax)
    If QApo >= 1000 Then Print QApo / 1000; "Km" Else Print QApo; "m"
    If QApo <= 0 Then 2010
    Print "     QPer         =";: QPer = Per(N1) - R(Mmax)
    If QPer >= 1000 Then Print QPer / 1000; "Km" Else Print QPer; "m"
    If QPer <= 0 Then 2010
    2010 Print "     Dist. media  =";
    If SM(N1) < 1000 Then Print SM(N1); "m"
    If SM(N1) < AL Then Print SM(N1) / 1000; "Km" Else Print SM(N1) / AL; "al"
    2020 Locate , 6
    Print "Vt =";
    If vt(N1) >= 1000 Then Print vt(N1) / 1000; "Km/s " Else Print vt(N1); "m/s "
    If nc > 2 Then 2030
    Print "     Va ="; DA / it; "rd/s "
    Print "     Vr ="; DA * (pdist(1, 2) + Dist(1, 2)) ^ 2 / 2 / it; "m2/s "
    2030 A$ = InKey$: If A$ = "" Then 2030
    If Asc(A$) > 47 And Asc(A$) < 58 Then PCopy 1, 0: GoTo 2000
    If A$ = " " Then PCopy 1, 0: GoSub globali: I2 = 1: GoTo 210
    2040 Return
    '
    globali:
    Locate 2, 1: Color 3
    Print "Tempo =";
    tanni = 0: tmesi = 0: tgiorni = 0: tore = 0: tminuti = 0: tsecondi = 0
    tanni = Int(t / 31536000)
    tmesi = Int((t - (tanni * 31536000)) / 2592000)
    tgiorni = Int((t - (tanni * 31536000) - (tmesi * 2592000)) / 86400)
    tore = Int((t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400)) / 3600)
    tminuti = Int((t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400) - (tore * 3600)) / 60)
    tsecondi = Int(t - (tanni * 31536000) - (tmesi * 2592000) - (tgiorni * 86400) - (tore * 3600) - (tminuti * 60))
    Print tanni; "yr"; tmesi; "mo"; tgiorni; "d"; tore; "h"; tminuti; "m"; tsecondi; "s ";
    Print "="; t; "sec "
    Print "Num. corpi ="; Nt
    Color 14: Print " En. potenziale    = "; Ep; "Joule "
    Color 4: Print " En. cinetica      = "; Ec
    Color 2: Print " En. totale        = "; Etot
    Color 15: Print " En. pot. media    = "; Epmed; "Joule "
    Print " En. cin. media    = "; Ecmed
    Print " En. tot. media    = "; Etmed
    Color 9: Print " Entropia (Etot/Ecin) = "; Sen
    Color 15: Print " Entr. media          = "; Senmed;
    3000 A$ = InKey$: If A$ = "" Then 3000
    If Asc(A$) > 47 And Asc(A$) < 58 Then PCopy 1, 0: GoTo 2000
    Return
    '
    comandi:
    Color 6: Locate L, C: Print "   Comandi:"
    Locate , C: Print "[e] elenco comandi"
    Locate , C: Print "[p] parametri"
    Locate , C: Print "[q] dim. quadro"
    Locate , C: Print "[d] distanza/quota"
    Locate , C: Print "[v] velocita'"
    Locate , C: Print "[g] angolo direz."
    Locate , C: Print "[a] assetto"
    Locate , C: Print "[i] intervallo tempo"
    Locate , C: Print "[t] inverti tempo"
    Locate , C: Print "[r] randomizza"
    Locate , C: Print "[n] di nuovo"
    Locate , C: Print "[c] pulisci schermo"
    Locate , C: Print "[</>] aum./rid. dim. quadro"
    Locate , C: Print "[" + Chr$(27) + Chr$(18) + Chr$(26) + "] sposta fin. schermo"
    Locate , C: Print "[h] halt"
    Locate L + 1, C + 31: Print "[1-0] dati analitici (per max. 10 corpi)"
    Locate , C + 31: Print "[ ] dati globali"
    Locate , C + 31: Print "[SH+1] mostra/nascondi diagr. energia"
    Locate , C + 31: Print "[SH+2] mostra/nascondi diagr. entropia"
    Locate , C + 31: Print "[F1] vista panoramica"
    Locate , C + 31: Print "[F2] assonometria"
    Locate , C + 31: Print "[F3] proiezioni x-y e x-z"
    Locate , C + 31: Print "[x] ruota x in senso antiorario"
    Locate , C + 31: Print "[X] ruota x in senso orario"
    Locate , C + 31: Print "[y] ruota y in senso antiorario"
    Locate , C + 31: Print "[Y] ruota y in senso orario"
    Locate , C + 31: Print "[z] ruota z in senso antiorario"
    Locate , C + 31: Print "[Z] ruota z in senso orario"
    Locate , C + 31: Print "[ESC] fine";
    GoSub 4000
    Return
    '
    4000 A$ = InKey$: If A$ = "" Then 4000
    If A$ = Chr$(27) Then Screen 0: End
    Return
    '
    unit1:
    Locate 4, 1: Print "unita' di massa:"
    Print "1) Masse terrestri "
    Print "2) Masse solari "
    Print "3) Kg "
    5000 A$ = InKey$: If A$ = "" Then 5000
    If Asc(A$) < 49 Or Asc(A$) > 52 Then 5000
    If A$ = "1" Then
        um# = MT: um$ = "mt"
    ElseIf A$ = "2" Then um# = MS: um$ = "ms"
    Else um# = 1: um$ = "Kg"
    End If
    Locate 9, 1: Print "unita' di raggio:"
    Print "1) Raggi terrestri "
    Print "2) Raggi solari "
    Print "3) Km "
    Print "4) Anni-luce "
    5010 A$ = InKey$: If A$ = "" Then 5010
    If Asc(A$) < 49 Or Asc(A$) > 52 Then 5010
    If A$ = "1" Then
        ur# = RT: ur$ = "rt"
    ElseIf A$ = "2" Then ur# = RS: ur$ = "rs"
    ElseIf A$ = "3" Then ur# = 1000: ur$ = "Km"
    Else ur# = AL: ur$ = "al"
    End If: PCopy 1, 0
    Return
    '
    unit2:
    PCopy 1, 0
    Locate 3, 1: Print "unita' di distanza:"
    Print "1) Raggi terrestri "
    Print "2) Raggi solari "
    Print "3) Km "
    Print "4) Anni-luce "
    5020 A$ = InKey$: If A$ = "" Then 5020
    If Asc(A$) < 49 Or Asc(A$) > 52 Then 5020
    If A$ = "1" Then
        ud# = RT: ud$ = "rt"
    ElseIf A$ = "2" Then ud# = RS: ud$ = "rs"
    ElseIf A$ = "3" Then ud# = 1000: ud$ = "Km"
    Else ud# = AL: ud$ = "al"
    End If: PCopy 1, 0
    Return
    '
    assetto:
    6000 Cls: Locate 2, 1
    Print " Precisazione distanza dal I astro "
    Print " (piano x-y): "
    Locate 5, 1
    gx = lo / dmax0(Dmax) / 2 / 2: Rg(Mmax) = R0(Mmax) * gx
    xx1(Mmax) = 320: yy1(Mmax) = 240: yy2(Mmax) = 240
    Circle (xx1(Mmax), yy1(Mmax)), Rg(Mmax), 15, , , 1
    For N1 = 1 To nc: If N1 = Mmax Then 6020
        Rg(N1) = R0(N1) * gx
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6010
        If yy1(N1) < 0 Or yy1(N1) > lv Then 6010
        Circle (xx1(N1), yy1(N1)), Rg(N1), 0, , , 1
        If N1 <> Mmax Then AZ(N1) = pigreco / 2: EL(N1) = 0
        6010 x0(N1) = dmax0(N1) * Cos(AZ(N1)) * Cos(EL(N1))
        If Abs(x0(N1)) < .1# Then x0(N1) = 0
        y0(N1) = dmax0(N1) * Sin(AZ(N1)) * Cos(EL(N1))
        If Abs(y0(N1)) < .1# Then y0(N1) = 0
        xx1(N1) = x0(N1) * gx + 320.5: yy1(N1) = y0(N1) * gx + 240.5
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6020
        If yy1(N1) < 0 Or yy1(N1) > lv Then 6020
        Circle (xx1(N1), yy1(N1)), Rg(N1), 15, , , 1
        Print "Dist."; N1; "="; dmax0(N1) / ud#(N1); ud$(N1);
        Print Space$(25 - Len(Str$(dmax0(N1) / ud#(N1))))
    6020 Next N1
    Locate 20, 1: Print " Confermi (s) ? "
    6030 A$ = InKey$: If A$ = "" Then 6030
    If A$ <> "n" Then 6060
    Dmax = 1: For N1 = 1 To nc: If N1 = Mmax Then 6050
        6040 Locate 4, 1: L = 4: C = 11
        Print "Dist."; N1; "="; dmax0(N1) / ud#(N1); ud$(N1);
        Print Space$(25 - Len(Str$(dmax0(N1) / ud#(N1))))
        GoSub scrivi: If N > 0 Then dmax0(N1) = N * ud#(N1) Else GoTo 6040
        If dmax0(N1) > dmax0(Dmax) Then Dmax = N1
    6050 Next N1: GoTo 6000
    6060 Cls: Locate 2, 1
    Print " Inserimento azimut (piano x-y): "
    Print " [" + Chr$(27) + Chr$(26) + "] - poi: [Invio]"
    Circle (xx1(Mmax), yy1(Mmax)), Rg(Mmax), 15, , , 1
    For N1 = 1 To nc: If N1 = Mmax Then 6130
        6070 If AZ(N1) < 0 Then AZ(N1) = AZ(N1) + 2 * pigreco
        If AZ(N1) > 2 * pigreco Then AZ(N1) = AZ(N1) - 2 * pigreco
        x0(N1) = dmax0(N1) * Cos(AZ(N1)) '* COS(EL(N1))
        If Abs(x0(N1)) < .1# Then x0(N1) = 0
        y0(N1) = dmax0(N1) * Sin(AZ(N1)) '* COS(EL(N1))
        If Abs(y0(N1)) < .1# Then y0(N1) = 0
        xx1(N1) = x0(N1) * gx + 320.5: yy1(N1) = y0(N1) * gx + 240.5
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6080
        If yy1(N1) < 0 Or yy1(N1) > lv Then 6080
        Circle (xx1(N1), yy1(N1)), Rg(N1), 15, , , 1
        6080 Locate 10, 1
        Print " Az."; N1; "= "; CInt(AZ(N1) * 180 / pigreco); " gradi     "
        A$ = InKey$
        If A$ = Chr$(13) Then 6090
        If Len(A$) < 2 Then 6070
        AK = Asc(Right$(A$, 1))
        If AK = 75 Then AZ(N1) = AZ(N1) - (AZ(N1) >= -2 * pigreco) * pigreco / 180
        If AK = 77 Then AZ(N1) = AZ(N1) + (AZ(N1) <= 2 * pigreco) * pigreco / 180
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6070
        If yy1(N1) < 0 Or yy1(N1) > lv Then 6070
        Circle (xx1(N1), yy1(N1)), Rg(N1), 0, , , 1
        GoTo 6070
        6090 Locate 9, 1: Print " Precisa: ": L = 10: C = 11
        GoSub scrivi: AZ(N1) = N * pigreco / 180
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6100
        If yy1(N1) < 0 Or yy1(N1) > lv Then 6100
        Circle (xx1(N1), yy1(N1)), Rg(N1), 0, , , 1
        6100 x0(N1) = dmax0(N1) * Cos(AZ(N1)) * Cos(EL(N1))
        If Abs(x0(N1)) < .1# Then x0(N1) = 0
        y0(N1) = dmax0(N1) * Sin(AZ(N1)) * Cos(EL(N1))
        If Abs(y0(N1)) < .1# Then y0(N1) = 0
        xx1(N1) = x0(N1) * gx + 320.5: yy1(N1) = y0(N1) * gx + 240.5
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6110
        If yy1(N1) < 0 Or yy1(N1) > lv Then 6110
        Circle (xx1(N1), yy1(N1)), Rg(N1), 15, , , 1
        6110 Locate 20, 1: Print " Confermi (s) ? "
        6120 A$ = InKey$: If A$ = "" Then 6120
        If A$ = "n" Then 6070
    6130 Next N1
    6140 Cls: Locate 2, 1
    Print " Inserimento elevazione (piano x-z): "
    Print " [" + Chr$(27) + Chr$(26) + "] - poi: [Invio]"
    Circle (xx1(Mmax), yy2(Mmax)), Rg(Mmax), 15, , , 1
    For N1 = 1 To nc: If N1 = Mmax Then 6210
        6150 If EL(N1) < 0 Then EL(N1) = EL(N1) + 2 * pigreco
        If EL(N1) > 2 * pigreco Then EL(N1) = EL(N1) - 2 * pigreco
        x0(N1) = dmax0(N1) * Cos(EL(N1))
        If Abs(x0(N1)) < .1# Then x0(N1) = 0
        z0(N1) = dmax0(N1) * Sin(EL(N1))
        If Abs(z0(N1)) < .1# Then z0(N1) = 0
        xx1(N1) = x0(N1) * gx + 320.5: yy2(N1) = z0(N1) * gx + 240.5
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6160
        If yy2(N1) < 0 Or yy2(N1) > lv Then 6160
        Circle (xx1(N1), yy2(N1)), Rg(N1), 15, , , 1
        6160 Locate 10, 1
        Print " El."; N1; "= "; CInt(EL(N1) * 180 / pigreco); " gradi     "
        A$ = InKey$
        If A$ = Chr$(13) Then 6170
        If Len(A$) < 2 Then 6150
        AK = Asc(Right$(A$, 1))
        If AK = 75 Then EL(N1) = EL(N1) - (EL(N1) >= -2 * pigreco) * pigreco / 180
        If AK = 77 Then EL(N1) = EL(N1) + (EL(N1) <= 2 * pigreco) * pigreco / 180
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6150
        If yy2(N1) < 0 Or yy2(N1) > lv Then 6150
        Circle (xx1(N1), yy2(N1)), Rg(N1), 0, , , 1
        GoTo 6150
        6170 Locate 9, 1: Print " Precisa: ": L = 10: C = 11
        GoSub scrivi: EL(N1) = N * pigreco / 180
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6180
        If yy2(N1) < 0 Or yy2(N1) > lv Then 6180
        Circle (xx1(N1), yy2(N1)), Rg(N1), 0, , , 1
        6180 x0(N1) = dmax0(N1) * Cos(EL(N1)) '* COS(AZ(N1))
        If Abs(x0(N1)) < .1# Then x0(N1) = 0
        z0(N1) = dmax0(N1) * Sin(EL(N1))
        If Abs(z0(N1)) < .1# Then z0(N1) = 0
        xx1(N1) = x0(N1) * gx + 320.5: yy2(N1) = z0(N1) * gx + 240.5
        If xx1(N1) < 0 Or xx1(N1) > lo Then 6190
        If yy2(N1) < 0 Or yy2(N1) > lv Then 6190
        Circle (xx1(N1), yy2(N1)), Rg(N1), 15, , , 1
        6190 Locate 20, 1: Print " Confermi (s) ? "
        6200 A$ = InKey$: If A$ = "" Then 6200
        If A$ = "n" Then 6150
        x0(N1) = x0(N1) * Cos(AZ(N1)) * Cos(EL(N1))
        If Abs(x0(N1)) < .1# Then x0(N1) = 0
        y0(N1) = dmax0(N1) * Sin(AZ(N1))
        If Abs(y0(N1)) < .1# Then y0(N1) = 0
        z0(N1) = z0(N1) * Cos(AZ(N1)) * Cos(EL(N1))
        If Abs(z0(N1)) < .1# Then z0(N1) = 0
    6210 Next N1
    pkx = kx: GoSub pargraf: GoSub calcint
    For N1 = 1 To nc: xx1(N1) = -1: yy1(N1) = -1: Next N1
    If I2 Then
        Locate 20, 1
        Print " Puoi premere -g- per cambiare subito l'angolo di direzione "
        6220 A$ = InKey$: If A$ = "" Then 6220
        If A$ = "g" Then
            I2 = 0: For N1 = 1 To nc: Lc = 21: GoSub angdir: Next N1
            I1 = 1: GoTo 190
        End If
    End If
    If kx <> pkx Then Screen 9, , 1, 0: Cls
    PCopy 1, 0
    Return
    '
    dimquad:
    7000 Locate 14, 1: L = 15: C = 20
    Print "dim. quadro (dist. max.": Print "dal primo astro) ="; dq
    GoSub scrivi: dq = N
    If dq < 1 Then dq = 1: GoTo 7000
    GoSub pargraf: GoSub calcint
    Return
    '
    velorb:
    v(1) = Sqr(G * (M0(2) ^ 2) / (M0(1) + M0(2)) / dmax0(2))
    v(2) = Sqr(G * (M0(1) ^ 2) / (M0(1) + M0(2)) / dmax0(2))
    Return
    '
    veloc:
    L = 10: C = 20
    Locate L, 1
    Print "veloc."; N1; "(m/sec) ="; v(N1);
    Print Space$(23 - Len(Str$(v(N1))))
    GoSub scrivi: v(N1) = N
    Return
    '
    angdir:
    If I2 Then
        Locate 20, 1: Print " Puoi premere -a- per cambiare subito l'assetto "
    End If
    C = 31
    L = Lc: Locate L, 1: Print Space$(40): Print Space$(40)
    Locate L, 1
    Print "direz. piano x-y"; N1; "(gradi) = "; dirxy(N1); "  "
    GoSub scrivi: dirxy(N1) = N
    L = Lc + 1: Locate L, 1
    Print "direz. piano x-z"; N1; "(gradi) = "; dirz(N1); "  "
    GoSub scrivi: dirz(N1) = N
    7001 If I2 Then
        A$ = InKey$: If A$ = "" Then 7001
        If A$ = "a" Then I2 = 0: GoSub assetto: I1 = 1: GoTo 210
    End If
    Return
    '
    astromax:
    Mmax = 1: For N1 = 1 To nc
        If M(N1) > M(Mmax) Then Mmax = N1
    Next N1: Return
    '
    qudist:
    Locate 12, 1: Print "quota(1) o distanza(2) ? (2) ";
    8000 A$ = InKey$: If A$ = "" Then 8000
    If A$ <> "1" Then A$ = "2"
    If A$ = "2" Then GoSub distmax: Return
    GoSub unit2: For N1 = 1 To nc
        ud#(N1) = ud#: ud$(N1) = ud$: dmax0(N1) = ud#(N1)
    Next N1
    8010 q0 = 0: Locate 13, 1
    L = 13: C = 14: Print "quota ("; ud$(1); ") = "; q0 / ud#(1);
    GoSub scrivi: q0 = N * ud#(1)
    If q0 <= 0 Then q0 = (dmax0(2) - R(1) - R(2)) * 2: GoTo 8010
    dmax0(2) = R(1) + q0: dmax0(1) = 0
    Dmax = 2
    x0(1) = 0: y0(1) = 0
    x0(2) = dmax0(2) * Cos(EL(2)) * Cos(AZ(2))
    If Abs(x0(2)) < .1# Then x0(2) = 0
    y0(2) = dmax0(2) * Cos(EL(2)) * Sin(AZ(2))
    If Abs(y0(2)) < .1# Then y0(2) = 0
    z0(2) = dmax0(2) * Sin(EL(2))
    If Abs(z0(2)) < .1# Then z0(2) = 0
    GoSub pargraf: GoSub calcint
    Return
    '
    distmax:
    For N1 = 1 To nc: If N1 = Mmax Then 9010
        GoSub unit2
        ud#(N1) = ud#: ud$(N1) = ud$
        9000 L = 13: C = 28: Locate 13, 1
        Print "dist."; N1; "dal I astro ("; ud$(N1); ") ="; dmax0(N1) / ud#(N1);
        GoSub scrivi: dmax0(N1) = N * ud#(N1)
        If dmax0(N1) <= 0 Then dmax0(N1) = (R(N1) + R(Mmax)) * 2: GoTo 9000
    9010 Next N1
    Dmax = 1: For N1 = 1 To nc
        If dmax0(N1) > dmax0(Dmax) Then Dmax = N1
    Next N1
    For N1 = 1 To nc: If N1 = Mmax Then 9020
        x0(N1) = dmax0(N1) * Cos(EL(N1)) * Cos(AZ(N1))
        If Abs(x0(N1)) < .1# Then x0(N1) = 0
        y0(N1) = dmax0(N1) * Cos(EL(N1)) * Sin(AZ(N1))
        If Abs(y0(N1)) < .1# Then y0(N1) = 0
        z0(N1) = dmax0(N1) * Sin(EL(N1))
        If Abs(z0(N1)) < .1# Then z0(N1) = 0
    9020 Next N1: GoSub pargraf: GoSub calcint
    Return
    '
    pargraf:
    ' - Calcolo parametri grafici -
    kx = lo / dmax0(Dmax) / dq / 2
    For N1 = 1 To nc: Rg(N1) = R(N1) * kx: Next
    primo = 1
    Return
    '
    calcint:
    vmax1 = 0
    For N1 = 1 To nc: If M0(N1) = 0 Then 10000
        If v(N1) > vmax1 Then vmax1 = v(N1)
    10000 Next N1
    it = Int(dmax0(Dmax) / vmax1 / 200)
    Return
    '
    interv:
    11000 Locate 16, 1: L = 16: C = 20
    Print "intervallo (sec) ="; it: GoSub scrivi: it = N
    If it = 0 Then it = 1: GoTo 11000
    Return
    '
    randomizza:
    For N1 = 1 To nc
        um$(N1) = "Kg": ur$(N1) = "Km": ud$(N1) = "Km"
        M0(N1) = Int(Rnd * MS * 10) + MT
        M(N1) = M0(N1): Mtot = Mtot + M(N1)
        R0(N1) = Int(Rnd * RS * 10) + RT: R(N1) = R0(N1)
        dirxy(N1) = Int(Rnd * 360): dirz(N1) = Int(Rnd * 360)
    Next N1
    GoSub astromax
    For N1 = 1 To nc
        dmax0(N1) = Int(Rnd * RS * 1000) + RS * 10
        AZ(N1) = Rnd * 2 * pigreco + Rnd * 2 * pigreco / N1
        EL(N1) = Rnd * 2 * pigreco + Rnd * 2 * pigreco / N1
        x0(N1) = dmax0(N1) * Cos(AZ(N1)) * Cos(EL(N1))
        y0(N1) = dmax0(N1) * Sin(AZ(N1)) * Cos(EL(N1))
        z0(N1) = dmax0(N1) * Sin(EL(N1))
        v(N1) = Sqr(G * (M0(Mmax) ^ 2) / (M0(N1) + M0(Mmax)) / dmax0(N1))
    12000 Next N1
    Dmax = 1: For N1 = 1 To nc
        If dmax0(N1) > dmax0(Dmax) Then Dmax = N1
    Next N1
    GoSub pargraf: GoSub calcint
    GoTo 190
    Return
    '
    zprof1:
    For N1 = 1 To nc: If M(N1) = 0 Then 13000
        prof(N1) = N1
    13000 Next N1
    For N1 = 1 To nc: If M(N1) = 0 Then 13002
        For N2 = N1 + 1 To nc: If M(N2) = 0 Then 13001
            If zz1(N1) < zz1(N2) Then Swap prof(N2), prof(N1)
        13001 Next N2
    13002 Next N1
    Return
    
    zprof2:
    For N1 = 1 To nc: If M(N1) = 0 Then 14000
        prof(N1) = N1
    14000 Next N1
    For N1 = 1 To nc: If M(N1) = 0 Then 14002
        For N2 = N1 + 1 To nc: If M(N2) = 0 Then 14001
            If zz2(N1) < zz2(N2) Then Swap prof(N2), prof(N1)
        14001 Next N2
    14002 Next N1
    Return
    '
    scrivi:
    Locate L, C: Input "", W: N$ = ""
    For k = C To C + 30: N$ = N$ + Chr$(Screen(L, k)): Next
    N = Val(N$): Return
    '
    End
    
  • Re: Programma in QBasic difettoso

    Altri linguaggi di programmazione

    https://www.iprogrammatori.it/forum-programmazione/altri-linguaggi/

Devi accedere o registrarti per scrivere nel forum
21 risposte