Robitex's Blog

Ideas in the web

Archivi Categorie: Asymptote

Spettri sismici – Parte 2


Un codice più flessibile

Continuiamo il perfezionamento del codice che avevamo scritto in precedenza per il disegno degli spettri elastici di progetto secondo il capitolo 3 delle nuove norme sulle costruzioni (D.M. 14/01/2008), dipendente dai tre parametri base a_g/g, F_0, T_C^* e relativamente agli stati limite di operatività SLO, di danno SLD, di salvaguardiadella vita SLV e di collasso SLC.

In Asymptote, il software che useremo per svolgere i calcoli e plottare i grafici, esiste la possibilità di creare nuovi tipi di dato che incorporano dati e funzioni, ovvero l’invenzione della programmazione ad oggetti.

Questa possibilità è offerta da struct, il costrutto che assomiglia a quello del C++, che ci permetterà di rappresentare le Categorie di sottosuolo e le categorie topografiche del sito.

E cominciando da quest’ultima. Serve un oggetto che chiameremo TopoCat che memorizzi il valore del coefficiente topografico e che lo restituisca per mezzo di una funzione:

// oggetto per le categorie topografiche
struct TopoCat {
   private real st;
   void operator init(real st){
      this.st = st;
   }
   real St(){
      return this.st;
   }
}

Abbiamo dichiarato come privato il campo st così che solo le funzioni pubbliche ne costituiscano l’interfaccia, ed inserito lo speciale operatore init che definisce il cosidetto costruttore, con cui l’utente creerà nuove istanze dell’oggetto. In questo caso il costruttore inizializza il campo interno accessibile con il riferimento this in dot notation.

Per instanziare un oggetto per ciascuna categoria topografica scriveremo così, creando i quattro tipi classici previsti dalla tabella 3.2.VI della norma:

// categorie topografiche classiche
TopoCat T1 = TopoCat(1.0);
TopoCat T2 = TopoCat(1.2);
TopoCat T3 = TopoCat(1.2);
TopoCat T4 = TopoCat(1.4);

Anche per le categorie di sottosuolo creiamo un oggetto apposito che viene instanziato con i coefficienti della tabella 3.2.V per fornire, sottoforma di funzioni naturalmente, i valori di Ss e Cc.

Da notare come il codice gestisce l’espressione di Ss: invece di riferirsi a campi interni dell’oggetto, i valori di ag/g e F0 devono semplicemente essere passati alle funzioni pubbliche (dette anche ‘metodi’ nel gergo della Object Oriented Programming).
Credo anche che fare la stessa cosa con un foglio di calcolo non sia così immediato.

// Categoria di sottosuolo
struct SoilCat {
   private real minSs;
   private real maxSs;
   private real k1;
   private real k2;
   private real h1;
   private real h2;

   // constructor
   void operator init(real mn, real mx, real k1, real k2, real h1, real h2){
      this.minSs = mn;
      this.maxSs = mx;
      this.k1 = k1;
      this.k2 = k2;
      this.h1 = h1;
      this.h2 = h2;
   }

   real Ss(real agg, real f0){// return the Ss value
      private real s;
      s = this.k1 - this.k2 * f0 * agg;
      if ( s < this.minSs ) {
         s = this.minSs;
      }
      if ( s > this.maxSs ) {
         s = this.maxSs;
      }
      return s;
   }
   
   real Cc(real tc){// return the Cc value
      return this.h1 * tc^this.h2;
   }
}

//
// istanziamento oggetti categorie sottosuolo (Tab. 3.2.V)
//
SoilCat A = SoilCat( 1.00 , 1.00 , 1.00 ,  0.00 , 1.00 ,  0.00);
SoilCat B = SoilCat( 1.00 , 1.20 , 1.40 , -0.40 , 1.10 , -0.20);
SoilCat C = SoilCat( 1.00 , 1.50 , 1.70 , -0.60 , 1.05 , -0.33);
SoilCat D = SoilCat( 0.90 , 1.80 , 2.40 , -1.50 , 1.25 , -0.50);
SoilCat E = SoilCat( 1.00 , 1.60 , 2.00 , -1.10 , 1.15 , -0.40);

Bene, notate come ovunque nel codice ci siano le dichiarazioni di tipo (il tipo corrispondente viene premesso al nome di ogni variabile), come richiede il linguaggio fortemente tipizzato di Asymptote.

Il codice finora scritto consente di rappresentare due diverse caratteristiche del sito. E’ naturale inglobare le due categorie in un ulteriore nuovo oggetto per rappresentare unitariamente il sito ed esporre i metodi per il disegno degli spettri.

Chiamiamo questo oggetto Spettro:

// Struct contenenti i parametri sismici del
// sito per lo Stati Limite considerato
// e funzioni grafiche e numeriche
//
// agrel = ag/g  ,  f0= F0  , tcstar = TC*
//
// q = fattore di struttura
//
struct Spettro {
   private real agrel;
   private real f0;
   private real tcstar;
  
   private real q;
   
   private real eta;
   
   private real smax;
   
   // periodi d'intervallo
   private real tb;
   private real tc;
   private real td;
   
   // ecco quì le caratteristiche del sito
   // rappresentate come campi privati
   private SoilCat soilC;
   private TopoCat topoC;
   
  // array dei punti notevoli
   pair[] keydot;
   
   // funzioni di definizione per intervalli dello spettro
   real s1(real T) {
      return  this.smax * ((T/this.tb) + (1 - T/this.tb)/(this.eta * this.f0));
    }

   real s3(real T){
      return this.smax * this.tc / T;
   }

   real s4(real T){
      return this.smax * (this.tc * this.td ) / T^2;
   }
   
  // costruttore
   void operator init(SoilCat s, TopoCat t, real q, real ag, real f0, real tc){
      this.soilC = s;
      this.topoC = t;
      this.q = q;
      this.agrel = ag;
      this.f0 = f0;
      this.tcstar = tc;

      this.eta = 1/q;

      this.tc = s.Cc(tc)*tc;
      this.tb = this.tc/3;
      this.td = 4.0 * ag + 1.60;

      this.smax = ag * s.Ss(ag,f0) * t.St() * (1/q) * f0;

      this.keydot = new pair[]{
        (      0 , this.s1(0)),
        (this.tb , this.smax),
        (this.tc , this.smax),
        (this.td , this.s4(this.td))
      };
   }
   guide plot(real xmax){
     return graph(this.s1 , 0 , this.tb)                   &
            (this.tb , this.smax) -- (this.tc , this.smax) &
            graph(this.s3 , this.tc , this.td)             &
            graph(this.s4 , this.td , xmax);
   }
}

Un commento anche per questo ultimo codice: la struct inizia con la dichiarazione dei campi privati e prosegue con quella delle funzioni matematiche dei rami dello spettro sismico ricavate dalla norma. La loro definizione deve precedere l’operatore init, poiché il costruttore ne fa uso per inizializzare le coordinate dei punti estremi dei rami dello spettro ai periodi Tb, Tc, e Td, salvandoli in un array chiamato keydot. Infine, il metodo plot che restituisce l’intero “percorso” del grafico.

L’oggetto non memorizza ulteriori metodi per facilitare il disegno dei grafici perché così facendo si lascia totale libertà di configurazione dei dettagli. In fondo la curva dello spettro è un oggetto indipendente dalla “tela” sulla quale verrà plottata.

Il codice “utente”

Ecco come può essere utilizzato il codice appena scritto. Consideriamo un sito su un terreno tipo B in pendio (T2). I quattro spettri elastici (q=1) con i vari parametri ag/g, F0 e Tc* saranno per esempio:

//
// un "oggetto struct" per singolo spettro
// istanziati con i dati sismici del sito
//
Spettro slo = Spettro(B, T2, 1.00, 0.04873, 2.462, 0.239);
Spettro sld = Spettro(B, T2, 1.00, 0.06107, 2.484, 0.251);
Spettro slv = Spettro(B, T2, 1.00, 0.15182, 2.402, 0.290);
Spettro slc = Spettro(B, T2, 1.00, 0.19475, 2.379, 0.298);

Il risultato

Non rimane che eseguire i plottaggi con il seguente codice, per ottenere in formato PDF il grafico completo pronto per essere inserito in relazione:

//
// plot degli spettri
//
size(16cm, 10cm, IgnoreAspect);
xlimits( 0 , 4.32);
ylimits( 0 , 0.72);

xaxis("$T$ (s)" ,
      BottomTop(),
      LeftTicks( Label(fontsize(8)) ,
               ptick=lightgray ,
               pTick=gray ,
               extend = true
               )
      );

yaxis("$S_{e}/g$" ,
      LeftRight(),
      RightTicks( Label(fontsize(8)),
                  pTick=lightgray,
                  extend = true
               )
      );

// plottaggio spettri (l'ordine è dall'alto verso il basso
// per la corrispondenza con la legenda
draw(slc.plot(4.2), red    + 1.80, "SLC");
draw(slv.plot(4.2), blue   + 1.80, "SLV");
draw(sld.plot(4.2), orange + 1.80, "SLD");
draw(slo.plot(4.2), green  + 1.80, "SLO");

// disegno dei punti notevoli
pen punto = linewidth(1.28bp);
dot(slc.keydot, punto);
dot(slv.keydot, punto);
dot(sld.keydot, punto);
dot(slo.keydot, punto);

// legenda
add(legend(), (3.1 , 0.60) , UnFill );
Spettri Sismici, fare click per scaricare il PDF

Spettri Sismici, fare click per scaricare il PDF

Cosa ci aspetta il futuro

Potremo aggiungere al grafico alcune informazioni di input per dare un tocco di professionalità, e pensare ad un modulo che calcoli i parametri ag/g, F0, Tc* del sito inserendone le coordinate geografiche. Chissà staremo a vedere.
Ciao

Spettri sismici Parte 1


Spettri sismici

Seismic Spectrum

Seismic Spectrum

La nuova normativa D.M. 14/01/2008 “Norme Tecniche per le Costruzioni”, in breve NTC 2008, entrata in vigore a metà del 2009 poco dopo il sisma in Abruzzo, ha introdotto una nuova mappatura del rischio sismico basata sul progetto S1 dell’INGV – Istituto Nazionale di Geofisica e Vulcanologia.

La pericolosità sismica dipende ora non più dai confini amministrativi dei comuni, ma dal punto d’interesse stesso, individuato in coordinate geografiche nel datum ED50.

Se anche a voi piace determinare gli spettri per proprio conto vi propongo un percorso tra gli strumenti del mondo (La)TeX, in particolare con Asymptote, con cui ottenere il grafico degli spettri sismici in formato vettoriale, pronto da inserire nella relazione tecnica del progetto strutturale.

Asymptote e gli spettri sismici

Oltre che ad imparare la programmazione con Asymptote (è consigliabile consultare il manuale di Asymptote presente nella vostra installazione o reperibile in rete), costruiremo gli spettri passo passo leggendo direttamente il paragrafo 3.2.3 della NTC 2008.

Questo post introduce l’argomento concentrandosi solo sullo spettro elastico di progetto allo stato limite di salvaguardia della vita per la sola componente orizzontale, e per un terreno di pianura di tipo B. Nei prossimi vedremo come il linguaggio di Asymptote possa implementare una pratica e completa interfaccia ai dati ed alle elaborazioni.

Installazione di Asymptote

Se avete una distribuzione LaTeX recente già installata sul vostro PC, qualsiasi sia il sistema operativo, tutto il software necessario è già pronto per l’uso, altrimenti potete rivolgervi indifferentemente a TeX Live 2009 (per Linux, Mac OS X, Windows, ecc) o MiKTeX 2.8 (per Windows), oppure potete installare esclusivamente Asymptote scaricandolo da questo link ed inserire nella relazione il grafico in uno dei nel formati grafici per esempio gif o png.

Calcolare lo spettro

Cominciamo!!! Avviate un editor di testo (Blocco Note, SciTE, Emacs, Gedit, ecc ecc), e scriviamo il seguente codice per instanziare le prime variabili:

// loading modulo graph
import graph;

// accelerazione di gravità (valore standard) in m/s2
real g = 9.80665;

// coefficiente di amplificazione topografica
real st = 1.0;

// fattore di struttura (spettri elastici)
real q = 1.0;

// eta per spettri di progetto
real eta = 1/q;

La sintassi è simile a Java, con le dichiarazione di tipo (real per le variabili numeriche in virgola mobile), il punto e virgola finale e le righe di commento che iniziano con il doppio slash.

Vediamo ora la prima funzione che restituisce il valore del coefficiente stratigrafico per un suolo di tipo B dalla tabella 3.2.V della NTC:

\displaystyle 1{,}00 \leq 1{,}40 - 0{,}40 F_0 \frac{a_g}{g}\leq 1{,}20

e la sua implementazione in Asymptote:

// funzione del coefficiente di amplificazione
// stratigrafica per il suolo tipo B
real ss( real agsug, real f0 ) {
	real tmp = 1.40 - 0.40 * f0 * agsug;
	if ( tmp < 1.00 ) {
		tmp = 1.0;
	}

	if ( tmp > 1.20 ) {
		tmp = 1.20;
	}
	return tmp;
}

Scriviamo adesso il codice per determinare i periodi caratteristici dello spettro a partire dai parametri fondamentali del sito a_g/g, F_0, T_C^*, per un periodo di ritorno corrispondente allo SLV, e per un sito di esempio:

// parametri fondamentali del sito
// agrel=ag/g, f0=F0, tcstar=Tc*
//
real agrel  = 0.1946;
real f0     = 2.423;
real tcstar = 0.28;

real tc = 1.10 * tcstar^(-0.20) * tcstar;
real tb = tc / 3 ;
real td = 4.0 * agrel + 1.60;

L’espressione dei quattro tratti dello spettro sono le seguenti (3.2.3.2.1 NTC 2008):

\displaystyle S_1 = a_g S \eta F_0 \left(\frac{T}{T_B}+\frac{1}{\eta F_0}\left(1-\frac{1}{T_B}\right)\right)

\displaystyle S_2 = a_g S \eta F_0

\displaystyle S_3 = a_g S \eta F_0\frac{T_c}{T}

\displaystyle S_4 = a_g S \eta F_0 \frac{T_cT_D}{T^2}

assegnamo quindi altrettante funzioni che restituiscono il valore dello spettro rapportato all’accelerazione di gravità:

// costante moltiplicativa
real s0 = agrel * ss(agrel, f0) * st * f0;

// primo tratto dello spettro
real s1(real T) {
  return  s0 * ((T/tb) + 1/(eta*f0) * (1 - (T/tb)));
}

// secondo tratto
real s2(real T) {
   return s0;
}

// terzo tratto
real s3(real T) {
   return s0 * tc / T;
}

// quarto tratto
real s4(real T) {
   return s0 * tc * td / T^2;
}

Graficare lo spettro

Tutto è pronto per disegnare il grafico. Basta scrivere una funzione chiamata appunto plot, che restituisce un oggetto di tipo guide, ovvero una curva a sua volta restituita dalla funzione di libreria graph. Il simbolo & è l’operatore di concatenazione delle singole curve, mentre il tratto orizzontale dello spettro è stato ottenuto tracciando direttamente un segmento tra il punto (tb, s0) ed il punto (tc, s0), utilizzando il simbolo -- secondo una sintassi molto comune tra i pacchetti grafici in LaTeX. La funzione graph accetta tre argomenti: il primo è la funzione che determina la curva mentre il secondo ed il terzo, delimitano l’intervallo di tracciamento sull’asse x.
La nostra funzione plot, accetta un numero inteso come limite complessivo del grafico sull’asse x, in modo da poter controllare la “coda”:

guide plot(real xmax){
	return 	graph(s1, 0, tb) &
		(tb, s0) -- (tc, s0) &
		graph(s3, tc, td)    &
		graph(s4, td, xmax);
}

Prepariamo l’area del grafico con l’ausilio di una libreria apposita integrata in Asymptote chiamata graph (modulo spiegato in dettaglio nel manuale di Asymptote):

// dimensione del grafico
size(16cm, 10cm, IgnoreAspect);

xlimits( 0 , 4.32);
ylimits( 0 , 0.70);

// etichettatura assi e caratterizzazione griglia
xaxis("$T$ (s)",
	BottomTop(),
	LeftTicks( Label(fontsize(8)),
		ptick=lightgray,
		pTick=gray,
		extend = true
	)
		);
yaxis("$S_{e}/g$" ,
	LeftRight(),
	RightTicks( Label(fontsize(8)),
		pTick=lightgray,
		extend = true
		)
	);

Plottiamo lo spettro passando alla funzione draw la nostra curva ed aggiungiamo una legenda:

// plottaggio spettro
draw(plot(4.2) , blue + 1.80 , "SLV");

// aggiungi legenda
add(legend(), (2.8 , 0.6) , UnFill );

Da notare come il linguaggio di Asymptote gestisca le proprietà di spessore e colore del tratto, semplicemente sommando il nome del colore con il valore dello spessore. Forse non sarà così elegante e qualche programmatore Java storcerà il naso, ma è senza dubbio comodo.
Tutte le funzioni sono spiegate nel manuale di Asymptote, in particolare la funzione add() e la funzione legend() utilizzate nell’ultima parte del codice.

Confezioniamo il tutto

Se avete seguito fin qui il codice digitato nel vostro editor dovrebbe essere questo:

// loading modulo graph
import graph;

// accelerazione di gravità (valore standard) in m/s2
real g = 9.80665;

// coefficiente di amplificazione topografica
real st = 1.0;

// fattore di struttura (spettri elastici)
real q = 1.0;

// eta per spettri di progetto
real eta = 1/q;

// funzione del coefficiente di amplificazione
// stratigrafica per il suolo tipo B
real ss( real agsug, real f0 ) {
	real tmp = 1.40 - 0.40 * f0 * agsug;
	if ( tmp < 1.00 ) {
		tmp = 1.0;
	}

	if ( tmp > 1.20 ) {
		tmp = 1.20;
	}
	return tmp;
}

//
// parametri fondamentali del sito
// agrel=ag/g, f0=F0, tcstar=Tc*
//
real agrel  = 0.1946;
real f0     = 2.423;
real tcstar = 0.28;

real tc = 1.10 * tcstar^(-0.20) * tcstar;
real tb = tc / 3 ;
real td = 4.0 * agrel + 1.60;

// costante moltiplicativa
real s0 = agrel * ss(agrel, f0) * st * f0;

// primo tratto dello spettro
real s1(real T) {
  return  s0 * ((T/tb) + 1/(eta*f0) * (1 - (T/tb)));
}

// secondo tratto
real s2(real T) {
   return s0;
}

// terzo tratto
real s3(real T) {
   return s0 * tc / T;
}

// quarto tratto
real s4(real T) {
   return s0 * tc * td / T^2;
}

guide plot(real xmax){
	return 	graph(s1, 0, tb)     &
		(tb, s0) -- (tc, s0) &
		graph(s3, tc, td)    &
		graph(s4, td, xmax);
}

// dimensione del grafico
size(16cm, 10cm, IgnoreAspect);

xlimits( 0 , 4.32);
ylimits( 0 , 0.70);

// etichettatura assi e caratterizzazione griglia
xaxis("$T$ (s)",
	BottomTop(),
	LeftTicks( Label(fontsize(8)),
		ptick=lightgray,
		pTick=gray,
		extend = true
	)
		);
yaxis("$S_{e}/g$" ,
	LeftRight(),
	RightTicks( Label(fontsize(8)),
		pTick=lightgray,
		extend = true
		)
	);

// plottaggio spettro
draw(plot(4.2) , blue + 1.80 , "SLV");

// aggiungi legenda
add(legend(), (2.8 , 0.6) , UnFill );

Bene. Salvato il file con l’estensione *.asy per esempio con il nome spettro.asy, basta aprire una finestra di comando (in gergo *nix il Terminale) e scrivere il comando di compilazione (spostatevi con il comando cd nella directory dove risiede il file):

asy -f pdf spettro.asy

Se non avviene automaticamente, aprite il file pdf di output e controllate il risultato. Confrontatelo anche con questo compilato da me.

Per oggi chiudiamo qui. La prossima volta ci occuperemo di generalizzare per le altre situazioni il codice con il paradigma di programmazione object oriented, concettualizzando lo spettro in un oggetto.
Alla prossima.

Il triangolo di Sierpiński


Il disegno della ricorsione

Visualizzare la ricorsione è l’idea di questo post. E la figura ricorsiva più semplice da disegnare mi sembra proprio il triangolo di Sierpiński. Dopo qualche notizia essenziale sul software che userò, passeremo al codice…

Sierpinski Triangle

Sierpinski Triangle of 6 level (click to download the pdf file)

Il linguaggio grafico di Asymptote

Per le nostre sperimentazioni useremo un software chiamato Asymptote, che consiste in un linguaggio di alto livello derivato dal METAPOST, piuttosto espressivo, simile a C++ e Java, quindi fortemente tipizzato.

Si possono instanziare i classici tipi numerici come int o real, e produrne di nuovi per mezzo del costrutto struct.

Infine, occorre terminare le istruzioni con il punto e virgola, utilizzare le parentesi graffe per definire i blocchi di codice, e produrre i commenti nel formato C++, ma vediamo subito un esempio (supporremo che una versione di Asymptote sia installata sul vostro sistema Windows, Linux o Mac OS che sia, magari tramite la TeX Live 2009).

Per disegnare una spezzata scrivete in un file di testo con estensione .asy il seguente codice, dove le coordinate sono separate da un doppio trattino che significa come in molti altri linguaggi del mondo LaTeX, una linea retta:

draw((0,0)--(2.5,2.5)--(5,0.5)--(7.5,5));

Poi in una finestra di console, date il comando di compilazione: asy -f pdf nomefile.asy. Al termine dell’elaborazione verrà rilasciato un file pdf contenente il risultato grafico desiderato.

Qualcosa di più complesso

Cominciamo a disegnare quello che sarà il triangolo di livello 0, utilizzando il tipo pair ovvero il tipo che in Asymptote corrisponde ad una coppia di numeri reali:

// triangle line length
real ell = 100;

// level 0 triangle:
pair p1 = (0,0);
pair p2 = (ell,0);
pair p3 = (ell/2,sqrt(3)*ell/2);
draw(p1--p2--p3--cycle);

Il termine cycle specifica che il percorso da disegnare sarà chiuso e corrisponde di fatto al primo punto p1.

Sierpiński

Il triangolo del matematico polacco venne da lui elaborato nel 1915, quando i sovietici tentavano di impossessarsi del cuore di quel popolo.

Disegnato un triangolo (che supporremo equilatero) detto di livello 0, si disegna il triangolo i cui vertici sono i punti di mezzo dei lati del triangolo principale, poi si applica ricorsivamente lo stesso procedimento ai tre sotto triangoli nord, sud-est e sud-ovest.

La struttura è un frattale in quanto è autosomigliante. Una semplice osservazione ci dice che man mano che si scende di livello le lunghezze dei lati dei triangoli si dimezzano, per cui se l_0 è la misura del lato del triangolo del livello 0, quella del lato a livello n sarà:

\displaystyle l_n = \frac{l_0}{2^n}.

Il numero dei triangoli t in funzione del livello v si ottiene considerando che scendendo di un livello il loro numero cresce di 3^v-1. In formule:

\displaystyle t = 1 + \sum_{i=1}^v 3^{i-1}

Per esempio al livello 5 i triangoli disegnati sono 122, ma al 10° salgono a 29.525 ed al 14° sono ben 2.391.485.

Il codice (non ottimizzato)

Per terminare i cicli ricorsivi altrimenti infiniti, è possibile utilizzare un criterio grafico: se lo spessore di tracciamento delle linee è tale che i triangoli dei livelli ulteriori risultano coperti, è inutile proseguire con il disegno.

La funzione ricorsiva chiamata sierpinski() prende i tre punti vertici del triangolo ed effettua immediatamente il test seguente (s è lo spessore del tratto ed l la lunghezza del lato del triangolo):

\displaystyle \frac{s}{2}\geq \frac{l}{2\sqrt{3}}

Se il test è falso significa che possiamo procedere al disegno del sotto-triangolo con i vertici nei punti di mezzo dei lati del triangolo principale, ed incrementare l’intero che conterrà alla fine il numero totale dei triangoli disegnati.

Il passo finale è lanciare ricorsivamente tre volte la funzione stessa sui sotto-triangoli del livello successivo.

Il codice non è ottimizzato perché in realtà non serve calcolare la distanza tra due punti per ottenere la lunghezza del lato del triangolo del livello attuale perché è ottenibile semplicemente dimezzando quella corrispondente al livello precedente.

Ecco il listato:

//
// Sierpinsky triangle: draw it with a recursive function
// Copyright (c) 2009 Roberto Giacomelli
// Term of licence see: https://robitex.wordpress.com/legalese/
//

// var setup
int n = 1;
real l0 = 100;
real s = 1.0;

// set the linewidth
defaultpen(s);

// main triangle
pair p1 = (0,0);
pair p2 = (l0,0);
pair p3 = (l0/2,sqrt(3)*l0/2);
draw(p1--p2--p3--cycle);

// function med() return the midpoint
// of the line ab
pair med(pair a, pair b){
return ((a.x+b.x)/2 , (a.y+b.y)/2);
}

// function dist() computes the length
// of line ab
real dist(pair a, pair b){
return sqrt((a.x-b.x)^2+(a.y-b.y)^2);
}

// main recursive function
void sierpinski(pair a, pair b, pair c){
// check graphic stop condition
if ( s<= dist(a,b)/sqrt(3) ){return;}

// draw subtriangle of the level
pair am = med(b,c);
pair bm = med(a,c);
pair cm = med(a,b);
draw(am--bm--cm--cycle);
++n;

//well, and now recursive running
// sublevel south east triangle
sierpinski(a,cm,bm);
// sublevel south west triangle
sierpinski(cm,b,am);
// sublevel north triangle
sierpinski(c,bm,am);
}
// draw the Sierpinski figure
sierpinski(p1,p2,p3);

// write on console the number
// of drawed triangles
write(n);

Scarica l’esempio in pdf : Sierpinski Triangle Figure
Alla prossima!

%d blogger hanno fatto clic su Mi Piace per questo: