StarDance - Come creare il programma

Come creare il programma

Il programma StarDance è stato realizzato utilizzando come struttura base un semplice programma N-corpi, il cui listato in QuickBasic è mostrato di seguito, che simula le interazioni gravitazionali esistenti in un gruppo di stelle. In più, StarDance utilizza un passo temporale "variabile", proporzionale alla minima distanza che intercorre fra due stelle qualsiasi del sistema, per affrontare il cosiddetto "problema degli incontri ravvicinati".

Il programma utilizza come unità di misura della distanza il parsec (1 pc = 3,09 * 10^18 cm), che corrisponde a 3,26 anni luce. La massa è invece misurata in masse solari (1 Ms = 1,99 * 10^33 g), mentre la misurazione del tempo è in milioni di anni (1 Ma = 3,16 * 10^13 s). Con questa scelta delle unità di misura la costante di gravitazione universale G assume il valore di 4,5 * 10^-3. E' comodo inoltre esprimere le velocità iniziali delle varie stelle in km/s; in tal caso, occorre moltiplicare i valori di queste velocità iniziali, espressi in km/s, per il fattore 1,02, che opera la conversione in unità pc/Ma.

REM **************************************************************************    
REM * * *  Programma N-corpi - Versione minima per PC IBM e compatibili  * * *
REM **************************************************************************
    DEFDBL A-Z
    SCREEN 12
REM -- Fisso i parametri della simulazione
    N = 2
    h = .5
    eps = 0
    tcrit = 200
    ingrand = 20
REM -- Inizializzo le matrici x, y, z, vx, vy, vz, m
    m(1) = 1
    m(2) = 1
    x(1) = 10
    y(1) = .5
    vx(1) = -.1 * 1.023
    x(2) = -10
    y(2) = -.5
    vx(2) = .1 * 1.023
REM -- Calcola l'evoluzione temporale del sistema
10  FOR i = 1 TO N
        ax(i) = 0
        ay(i) = 0
        az(i) = 0
        FOR j = 1 TO N
            IF j = i THEN 20
REM -- Calcola distanza e forza tra le stelle i e j
            d = SQR((x(i) - x(j)) ^ 2 + (y(i) - y(j)) ^ 2 + (z(i) - z(j)) ^ 2)
            f = (4.5 * 10 ^ -3) * m(j) / (d ^ 2 + eps ^ 2) ^ 1.5
REM -- Calcola l'accelerazione della stella i-esima
            ax(i) = ax(i) + f * (x(j) - x(i))
            ay(i) = ay(i) + f * (y(j) - y(i))
            az(i) = az(i) + f * (z(j) - z(i))
20      NEXT j
REM -- Aggiorna la velocità della stella i-esima
        vx1(i) = vx(i) + h * ax(i)
        vy1(i) = vy(i) + h * ay(i)
        vz1(i) = vz(i) + h * az(i)
REM -- Aggiorna la posizione della stella i-esima
        x1(i) = x(i) + h * vx1(i)
        y1(i) = y(i) + h * vy1(i)
        z1(i) = z(i) + h * vz1(i)
REM -- Visualizza le stelle sullo schermo
        PSET (x1(i) * ingrand + 320, 265 - y1(i) * ingrand), i + 3
    NEXT i
REM -- Assegna le nuove condizioni iniziali
    FOR i=1 TO N
        vx(i) = vx1(i)
        vy(i) = vy1(i)
        vz(i) = vz1(i)
        x(i) = x1(i)
        y(i) = y1(i)
        z(i) = z1(i)
    NEXT i
REM -- Calcola e visualizza il tempo fittizio
    tempo = tempo + h
    IF tempo > tcrit THEN GOTO 100
    LOCATE 3, 6: PRINT "TEMPO  = "; INT(tempo)
    LOCATE 3, 21: PRINT "milioni di anni "
    GOTO 10
100 END