Mandelbrot en pseudo 3D

'              3D_Mendelbrot.bas
dim x,y,z,u,v,w,a,b,l,sa,sb,xmax,xmin,ymax,ymin,fin,seuil
caption 0, " <CLICK> pour arrêter"
width 0,700 : height 0,500
xmax = 1.8 : xmin = -2.5 : ymax = 1.8 : ymin = -1.8 : fin = 65 : seuil = 4
2d_fill_color 0,0,255 : 2d_rectangle 0,0,1000,700
for v = 180 to 380        : '  Si vous êtes pressé ajouter step 2 ou toute autre valeur
  for u = 380-v to 820-v  : '  step 2
    y =(ymax-ymin)*(380-v)/200+ymin
    x =(xmax-xmin)*(u+v-380)/440+xmin
    z = 0 : a = 0 : b = 0 : l = 0
    while (l <= seuil)and(z <= fin)
      z = z + 1 : sa = a : sb = b : a = sa*sa-sb*sb+x : b = 2*sa*sb+y : l = a*a+b*b
    end_while
    w = z*3
    2d_pen_color 255,0,0    : 2d_line u,v,u,v-w
    2d_pen_color 0,0,0      : 2d_line u-2,v,u-2,v-w
    2d_pen_color 255,255,255 : 2d_line u-2,v-w,u,v-w
    if scancode <> 0 then terminate
  next u
next v
caption 0 ,"Terminé"

'              3D_Mendelbrot.bas
dim x,y,z,u,v,w,a,b,l,sa,sb,xmax,xmin,ymax,ymin,fin,seuil
caption 0, " <CLICK> pour arrêter"
width 0,720 : height 0,440
xmax = 1.8 : xmin = -2.5 : ymax = 1.8 : ymin = -1.8 : fin = 65 : seuil = 4
2d_fill_color 0,0,255 : 2d_rectangle 0,0,1000,700
for v = 180 to 380
  for u = 380-v to 820-v
    y =(ymax-ymin)*(380-v)/200+ymin
    x =(xmax-xmin)*(u+v-380)/440+xmin
    z = 0 : a = 0 : b = 0 : l = 0
    while (l <= seuil)and(z <= fin)
      z = z + 1 : sa = a : sb = b : a = sa*sa-sb*sb+x : b = 2*sa*sb+y : l = a*a+b*b
    end_while
    w = z*3
    2d_pen_color 30 + abs(440 -u)/4 ,0,0    : 2d_line u,v,u,v-w
    2d_pen_color abs(440 -u)/4,50,(v-180)/3      : 2d_line u-2,v,u-2,v-w
    2d_pen_color 30 + abs(440 -u)/2,155,255-(v-180)/3 : 2d_line u-2,v-w,u,v-w
    if scancode <> 0 then terminate
  next u
next v
caption 0 ,"Terminé"

rem ****************************************************************************
rem *
rem *
rem *      Adaptation d'un programme écrit en GWBASIC
rem *
rem *                    PAYSAGES FRACTALS
rem *
rem *                        PAR PAPYDALL
rem *
rem *
rem ****************************************************************************

label init,generer_matrice,calc_altitude,calc_centrage,calc_coord_ecran
label calc_coord_spatiales,affichage ,trace_contours ,remplis_facettes
label coord_ecran,test_phi,calc1,calc2,calc3,calc4,calc5,calc6,calc7,calc8,calc9
label suite1,suite2,coloriage,demarrer,fin
' ------------------------------------------------------------------------------
dim v1(3),v2(3),vn(3),vecl(3),vobs(3),pt(3,3),xe(3),ye(3),bary(3),c(3),type,r$
dim pi,coef,fang,mail,profil,mer,theta,phi,ray,ct,st,cp,sp,xobs,yobs,zobs
dim alfa,beta,r,xecl,yecl,zecl,dn,pas,ech,i,j,a,b,cc,d,alt,xemin,yemin
dim xemax,yemax,stp,ci,cj,xxe,yye,ecrx,ecry,rap,echx,echy,xcent,ycent
dim deb,fin,sens,l,tst,k,indm,ncoul,vvn,n,x1,x2,y1,prosc,vvobs
dim cosang,ang,freq,exx,exy,ex1,ex2,ex3,ey1,ey2,ey3
dim stp1,stp2,stp3,pch,comp,cpt1,cpt2,vvecl,coul,bord$
dim titre$,nl$ : nl$=chr$(13)
' ------------------------------------------------------------------------------

gosub demarrer : gosub fin
end

' ******************************************************************************
demarrer:
  gosub init : gosub generer_matrice : gosub calc_centrage : gosub affichage
return

 
' ******************************************************************************
init:
  Application_title "PAYSAGES FRACTALS"
  titre$ = " GENERATEUR DE PAYSAGES FRACTALS  "
  pi = 4 *atn(1) : coef = pi/180 : fang = pi/15
  width 0, 1000 : height 0,700
  color 0,50,100,255
 
' Le degré de maillage définit la taille des facettes élémentaires qui composent le paysage
' Plus ce degré est élevé et plus la représentation s'affine avec un réalisme de plus en plus grand
' mais le temps de calcul devient excessivement long.
  repeat
    r$ = message_input$("Saisie du degré de maillage","Degré de maillage (0 à 9)"+nl$+"mail = ","7")
  until  numeric(r$) = 1
  mail = val(r$)
' L'indice de profil correspond à la tendance qu'aura le relief à s'élancer vers le haut ou vers le bas
' Un grand indice de profil produira un relief très accidenté.
  repeat
    r$ = message_input$("Saisie de l'indice de profil","Indice de profil (1 à 100)"+nl$+"profil = ","100")
  until  numeric(r$) = 1
  profil = val(r$)
' Elévation du niveau de la mer : le choix de cette valeur détermine le moment où les parties dessinées
' seront considérées comme étant au-dessous du niveau de la mer.
  repeat
    r$ = message_input$("Saisie élévation de la mer","Elévation du niveau de la mer (-200 à 1000)"+nl$+"mer = ","500")
  until numeric(r$) = 1
  mer = val(r$)
' Angle de vision (vertical ou horizontal) et distance de l'observateur
' permettent de définir l'emplacement de l'observateur
  repeat
    r$ = message_input$("Saisie de l'angle de vision vertical","Angle de vision vertical theta (-90° à 90°)"+nl$+"theta = ","40")
  until numeric(r$) = 1
  theta = val(r$)
 
  repeat
    r$ = message_input$("Saisie de l'angle de vision horizontal","Angle de vision horizontal phi (0° à 360°)"+nl$+"phi = ","350")
  until numeric(r$) = 1
  phi = val(r$)
 
  repeat
    r$ = message_input$("Saisie distance","Distance de l'observateur à l'origine ( > 5000 )"+nl$+"ray = ","10000")
  until numeric(r$) = 1
  ray = val(r$)
 
  ct = cos(theta*coef) : st = sin(theta*coef)
  cp = cos(phi*coef)  : sp = sin(phi*coef)
  xobs = ray*ct*cp    : yobs = ray*ct*sp  : zobs = ray*st
 
' Types du tracé :
' 1- tracé rapide: le paysage apparaitra en tracé fil-de-fer,sans élimination des parties cachées
' 2- parties cachées : même type de tracé, mais avec élimination des parties cachées
' 3- Surfaces éclairées : le programme demande la position de la source lumineuse
  repeat
    r$ = message_input$("Saisie type de tracé","Type de tracé (1, 2, 3)"+nl$+"1. Rapide"+nl$+"2. Parties cachées"+nl$+"3. Surface éclairée"+nl$+"type = ","3")
  until numeric(r$) = 1
  type = val(r$)
 
  if type = 3
    repeat
      r$ = message_input$("Saisie de l'angle vertical d'éclairage","Angle vertical d'éclairage alfa (-90° à 90°)"+nl$+"alfa = ","50")
    until numeric(r$) = 1
    alfa = val(r$)*coef
    repeat
      r$ = message_input$("Saisie de l'angle horizontal d'éclairage","Angle horizontal d'éclairage beta (0° à 360°)"+nl$+"beta = ","250")
    until numeric(r$) = 1
    beta = val(r$)*coef
    repeat
      r$ = message_input$("Saisie de la source lumineuse","Distance de la source lumineuse à l'origine ( > 5000)"+nl$+"r = ","10000")
    until numeric(r$) = 1
    r = val(r$)
    xecl = r*cos(alfa)*cos(beta)
    yecl = r*cos(alfa)*sin(beta)
    zecl = r*sin(alfa)
    repeat
      r$ = message_input$("Saisie des contours des facettes","Désirez-vous les contours des facettes (O/N)"+nl$+"Bord$ = ","O")
    until upper$(r$) = "O" or upper$(r$) = "N"
    bord$ = r$
  end_if
  caption 0,titre$ + "  **** !!!  VEUILLEZ    PATIENTER  !!!  **** " + "  <ESC>  POUR ARRETER"
return
' ******************************************************************************
generer_matrice:
  dn = power(2,mail)+1 : pas = dn - 1 : ech = 4000
  dim noeud(dn,dn)
  while pas > 1
    for i = 1 to dn-pas step pas
      for j = 1 to dn-i-pas+1 step pas
        a = i+pas/2 : b = j+pas/2 : cc = i+pas : d = j+pas : gosub calc_altitude
        noeud(i,b) = (noeud(i,j)+noeud(i,d))/2+alt        : gosub calc_altitude
        noeud(a,j) = (noeud(i,j)+noeud(cc,j))/2+alt        : gosub calc_altitude
        noeud(a,b) = (noeud(cc,j)+noeud(i,d))/2+alt
      next j
    next i
    pas = pas/2 : ech = ech/2
  end_while
  return
' ******************************************************************************
' *** Calcul aléatoire des altitudes
calc_altitude:
  alt = rnd(1)*ech
  if rnd(1) > profil/100 then alt = 0-alt
return
' ******************************************************************************
' *** Calcul du centrage
calc_centrage:

  xemin = 1000 : yemin = 1000 : xemax = 0-1000 : yemax = 0-1000
  if mail > 5 then stp = power(2,(mail-5)) : else : stp = 1
  for i = 1 to dn step stp
    for j = 1 to dn-i+1 step stp
      ci = i : cj = j
      gosub calc_coord_ecran
      if xxe < xemin then xemin = xxe
      if xxe > xemax then xemax = xxe
      if yye < yemin then yemin = yye
      if yye > yemax then yemax = yye
    next j
  next i

  ecrx = 990 : ecry = 690  : ' Dimension de l'écran graphique
  rap = 1 : ' 2.62 : ' permet de tenir compte de la déformation d'écran
  echx = ecrx /(xemax-xemin)/rap
  echy = ecry /(yemax-yemin)
  if echy < echx then echx = echy
  echy = echx : echx = echx*rap
  xcent = (ecrx-echx*(xemax+xemin))/2
  ycent = (ecry+10-echy*(yemax+yemin))/2
return
' ******************************************************************************
affichage:
  if type = 1 then gosub calc3 : return
  if phi >= 120 and phi < 240 then gosub test_phi :return
  if phi < 120 then gosub calc5 : return
  gosub calc4
return
' ******************************************************************************
' *** Calcul des coordonnées d'écran
calc_coord_ecran:
  gosub calc_coord_spatiales
coord_ecran:
  d = c(1)*cp*ct+c(2)*sp*ct+c(3)*st-ray
  xxe = (c(1)*sp-c(2)*cp)/d : yye = (c(1)*cp*st+c(2)*sp*st-c(3)*ct)/d
  yye = 0-yye
return
' ******************************************************************************
' *** Calcul des coordonnées spatiales
calc_coord_spatiales:
  c(1) = ((1-ci)/(dn-1)+1/3)*4000*sqr(3)
  c(2) = ((cj-1)+(ci-1)/2)*8000/(dn-1)-4000
  c(3) = noeud(ci,cj)
  if c(3) < mer then c(3) = mer
return
' ******************************************************************************
test_phi:
  if phi < 180 then deb = 1 : sens = 1 : else : fin = 1 : sens = 0-1
  for i = 2 to dn
    tst = 0
    if phi < 180 then fin = dn-i+1: else : deb = dn-i+1
    for j = deb to fin step sens
      if phi < 180 then gosub calc1: gosub calc2 : else : gosub calc2: gosub calc1
    next j
  next i
return
' ******************************************************************************
calc1:
  ci = i-1 : cj = j : gosub calc_coord_spatiales
  for l = 1 to 3
    v1(l) = c(l): pt(1,l) = c(l)
  next l
  ci = i : gosub calc_coord_spatiales
  for l = 1 to 3
    v1(l) = v1(l)-c(l): pt(2,l) = c(l)
  next l
  if tst = 0 then tst = 1 : return
  if phi < 180 then cj = j-1 : else : cj = j+1 : ci = i-1
  gosub calc_coord_spatiales
  for l = 1 to 3
    pt(3,l) = c(l)
  next l
  gosub remplis_facettes
return
' ******************************************************************************
calc2:
ci = i-1 : cj = j+1 : gosub calc_coord_spatiales
for l = 1 to 3
  v2(l) = c(l) : pt(1,l) = c(l)
next l
ci = i : cj = j : gosub calc_coord_spatiales
for l = 1 to 3
  v2(l) = v2(l)-c(l) : pt(2,l) = c(l)
next l
if tst = 0 then tst = 1 : return
if phi < 180 then ci = i-1 : else : cj =j+1
gosub calc_coord_spatiales
for l = 1 to 3
  pt(3,l) = c(l)
next l
gosub remplis_facettes
return
' ******************************************************************************
calc3:
  cc = 1
  for  i = 2 to dn
    for j = 1 to dn-i+1
      ci = i : cj = cj : gosub calc_coord_ecran
      xe(1) = xxe*echx+xcent : ye(1) = yye*echy+ycent
      ci = i-1 : gosub calc_coord_ecran
      xe(2) = xxe*echx+xcent : ye(2) = yye*echy+ycent
      cj = j+1 : gosub calc_coord_ecran
      xe(3) = xxe*echx+xcent : ye(3) = yye*echy+ycent
      gosub trace_contours
    next j
  next i
return
' ******************************************************************************
calc4:
  if phi < 300 then deb = 1 : sens = 1 : else : fin = 1 : sens = 0-1
  for i = dn-1 to 1 step -1
    tst = 0 : if phi < 300 then fin = i : else : deb = i
    for j = deb to fin step sens
      k = i+1-j
      if phi < 300 then gosub calc6 : gosub calc7 : else : gosub calc7 : gosub calc6
    next j
  next i
return
' ******************************************************************************
' ******************************************************************************
calc5:
  if phi > 60 then deb = 1 : sens = 1 : else : fin = 1 : sens = 0-1
  for i = 2 to dn
    tst = 0 : if phi > 60 then fin = dn-i+1 : else : deb = dn-i+1
    for j = deb to fin step sens
      if phi > 60 then gosub calc9 : gosub calc8 : else : gosub calc8 : gosub calc9
    next j
  next i
return
' ******************************************************************************
calc6:
  ci = j : cj = k + 1 : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = c(l) : pt(1,l) = c(l) : next l
  cj = k : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = v1(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi < 300 then ci = j-1 : cj = k+1 : else : ci = j+1
  gosub calc_coord_spatiales
  for l = 1 to 3 : pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
calc7:
  ci = j+1 : cj = k : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = c(l) : pt(1,l) = c(l) : next l
  ci = j : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = v2(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi < 300 then cj = k+1 : else : ci = j+1 : cj = k-1
  gosub calc_coord_spatiales
  for l = 1 to 3 : pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
calc8:
  ci = j+1 : cj = i-1 : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = c(l) : pt(1,l) = c(l) : next l
  ci = j : cj = i : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = v1(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi > 60 then cj = i-1 : else : ci = j+1
  gosub calc_coord_spatiales
  for l = 1 to 3 :  pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
calc9:
  ci = j : cj = i-1 : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = c(l) : pt(1,l) = c(l) : next l
  cj = i : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = v2(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi > 60 then ci = j-1 : else : ci = j+1 : cj = i-1
  gosub calc_coord_spatiales
  for l = 1 to 3 : pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
' *** Tracé des contours
trace_contours:
  if cc = 1 then 2d_pen_color 0,0,0 : else : 2d_pen_color 0,255,0
  2d_line xe(1),ye(1),xe(2),ye(2) : ' 2d_poly_to xe(3),ye(3) : 2d_poly_to xe(1),ye(1)
  2d_line xe(2),ye(2),xe(3),ye(3)
  2d_line xe(3),ye(3),xe(1),ye(1)

return
' ******************************************************************************
' *** Calcul des couleurs et remplissagedes facettes
remplis_facettes:
  indm = 0 : ncoul = 0
  if pt(1,3) = mer and pt(2,3) = mer and pt(3,3) = mer then indm = 1 : goto suite1
  vn(1) = v1(2)*v2(3)-v1(3)*v2(2)
  vn(2) = v1(3)*v2(1)-v1(1)*v2(3)
  vn(3) = v1(1)*v2(2)-v1(2)*v2(1)
  vvn = sqr(vn(1)*vn(1)+vn(2)*vn(2)+vn(3)*vn(3))
  for n = 1 to 3 : bary(n) = (pt(1,n)+pt(2,n)+pt(3,n))/3 : next n
  vobs(1) = bary(1)-xobs : vobs(2) = bary(2)-yobs: vobs(3) = bary(3)-zobs
  prosc = vn(1)*vobs(1)+vn(2)*vobs(2)+vn(3)*vobs(3)
  vvobs = sqr(vobs(1)*vobs(1)+vobs(2)*vobs(2)+vobs(3)*vobs(3))
  cosang = prosc/(vvn*vvobs)
  ang = atn(sqr(1-cosang*cosang)/cosang)
  if ang > 0 then ncoul = 1
  if ncoul = 1 or type = 2  then goto suite1
  vecl(1) = bary(1)-xecl : vecl(2) = bary(2)-yecl : vecl(3) = bary(3)-zecl
  prosc = vn(1)*vecl(1)+vn(2)*vecl(2)+vn(3)*vecl(3)
  vvecl = sqr(vecl(1)*vecl(1)+vecl(2)*vecl(2)+vecl(3)*vecl(3))
  cosang = prosc/(vvn*vvecl)
  ang = atn(sqr(1-cosang*cosang)/cosang)
  if ang < 0 then ang = ang + pi
  coul = int(ang/fang+.5)
  if coul <> 0 then freq = 15/coul : else : freq = power(10,30)
suite1:
  for n = 1 to 3
    c(1) = pt(n,1) : c(2) = pt(n,2) : c(3) = pt(n,3)
    gosub coord_ecran
    xe(n) = int(xxe * echx + xcent + .5) : ye(n) = int(yye * echy + ycent + .5)
  next n
  if ye(2) >= ye(1) and ye(2) >= ye(3)
    exx = xe(1) : exy = ye(1) : xe(1) = xe(2) : ye(1) = ye(2)
    xe(2) = exx : ye(2) = exy : goto suite2
  end_if
  if ye(3) >= ye(1) and ye(3) >= ye(2)
    exx = xe(1) : exy = ye(1) : xe(1) = xe(3) : ye(1) = ye(3)
    xe(3) = exx : ye(3) = exy
  end_if
suite2:
  if ye(3) > ye(2)
    exx = xe(2) : exy = ye(2) : xe(2) = xe(3)
    ye(2) = ye(3) : xe(3) = exx : ye(3)= exy
  end_if
  ex1 = xe(1)-xe(2) : ey1 = ye(1)-ye(2)
  ex2 = xe(1)-xe(3) : ey2 = ye(1)-ye(3)
  ex3 = xe(2)-xe(3) : ey3 = ye(2)-ye(3)
  if ey1 <> 0 then stp1 = ex1/ey1
  if ey2 <> 0 then stp2 = ex2/ey2
  if ey3 <> 0 then stp3 = ex3/ey3
  pch = int(rnd(1)*freq+1.5) : comp = 1
  cpt2 = 0
  for cpt1 = 0 to ey1
    x1 = int(xe(1)-cpt1*stp1+.5) : x2 = int(xe(1)-cpt2*stp2+.5): y1 = ye(1)-cpt1
    if ey1 = 0 then x1 = xe(2)
    if ey2 = 0 then x2 = xe(3)
    gosub coloriage
  next cpt1
  for cpt1 = 1 to ey3
    x1 = int(xe(2)-cpt1*stp3+.5) : x2 = int(xe(1)-cpt2*stp2+.5): y1 = ye(1)-cpt2
    gosub coloriage
  next cpt1
  if (type=2 and ncoul=0 and indm=0)or(type=3 and ncoul = 1)or upper$(bord$)="O"
    cc = 1 : gosub trace_contours
  end_if
return

' ******************************************************************************
coloriage:
  if scancode = 27  then caption 0,"Arrêté par l'utilisateur" : end
  if (ncoul = 1 and type = 2) or (mod(y1,2) = 0 and indm = 1)
    cc = 1 : else : cc = 0
  end_if
  if type = 2 or indm = 1 or ncoul = 1
    if cc = 1 then 2d_pen_color 0,0,0 : else : 2d_pen_color 0,255,0
    2d_line x1,y1,x2,y1
    cpt2 = cpt2 +1 : return
  end_if
  if x1 = x2
    if comp = int(pch+.5)
      2d_pen_color 0,0,0 : else : 2d_pen_color 0,255,0
    end_if
    2d_point x1,y1
    cpt2 = cpt2 +1 : return
  end_if
  for n = x1 to x2 step sgn(x2-x1)
    cc = 0
    if comp = int(pch+.5) then cc = 1 : pch = pch + freq
    if cc = 1 then 2d_pen_color 0,0,0 : else : 2d_pen_color 0,255,0
    2d_point n,y1
    comp = comp + 1
  next n
  cpt2 = cpt2 + 1
return
' *************************************************************************
fin:
Caption 0,"Terminé"
end
' ******************************************************************************

rem ****************************************************************************
rem *
rem *
rem *      Adaptation d'un programme écrit en GWBASIC
rem *
rem *                    PAYSAGES FRACTALS
rem *
rem *                        PAR PAPYDALL
rem *
rem *
rem ****************************************************************************

label init,generer_matrice,calc_altitude,calc_centrage,calc_coord_ecran
label calc_coord_spatiales,affichage ,trace_contours ,remplis_facettes
label coord_ecran,test_phi,calc1,calc2,calc3,calc4,calc5,calc6,calc7,calc8,calc9
label suite1,suite2,coloriage,demarrer,fin
' ------------------------------------------------------------------------------
dim v1(3),v2(3),vn(3),vecl(3),vobs(3),pt(3,3),xe(3),ye(3),bary(3),c(3),type,r$
dim pi,coef,fang,mail,profil,mer,theta,phi,ray,ct,st,cp,sp,xobs,yobs,zobs
dim alfa,beta,r,xecl,yecl,zecl,dn,pas,ech,i,j,a,b,cc,d,alt,xemin,yemin
dim xemax,yemax,stp,ci,cj,xxe,yye,ecrx,ecry,rap,echx,echy,xcent,ycent
dim deb,fin,sens,l,tst,k,indm,ncoul,vvn,n,x1,x2,y1,prosc,vvobs
dim cosang,ang,freq,exx,exy,ex1,ex2,ex3,ey1,ey2,ey3
dim stp1,stp2,stp3,pch,comp,cpt1,cpt2,vvecl,coul,bord$
dim titre$,nl$ : nl$=chr$(13)
' ------------------------------------------------------------------------------

gosub demarrer : gosub fin
end

' ******************************************************************************
demarrer:
  gosub init : gosub generer_matrice : gosub calc_centrage : gosub affichage
return


' ******************************************************************************
init:
  Application_title "PAYSAGES FRACTALS"
  titre$ = " GENERATEUR DE PAYSAGES FRACTALS  "
  pi = 4 *atn(1) : coef = pi/180 : fang = pi/15
  width 0, 1000 : height 0,700
  color 0,50,100,255

' Le degré de maillage définit la taille des facettes élémentaires qui composent le paysage
' Plus ce degré est élevé et plus la représentation s'affine avec un réalisme de plus en plus grand
' mais le temps de calcul devient excessivement long.
  repeat
    r$ = message_input$("Saisie du degré de maillage","Degré de maillage (0 à 9)"+nl$+"mail = ","7")
  until  numeric(r$) = 1
  mail = val(r$)
' L'indice de profil correspond à la tendance qu'aura le relief à s'élancer vers le haut ou vers le bas
' Un grand indice de profil produira un relief très accidenté.
  repeat
    r$ = message_input$("Saisie de l'indice de profil","Indice de profil (1 à 100)"+nl$+"profil = ","100")
  until  numeric(r$) = 1
  profil = val(r$)
' Elévation du niveau de la mer : le choix de cette valeur détermine le moment où les parties dessinées
' seront considérées comme étant au-dessous du niveau de la mer.
  repeat
    r$ = message_input$("Saisie élévation de la mer","Elévation du niveau de la mer (-200 à 1000)"+nl$+"mer = ","500")
  until numeric(r$) = 1
  mer = val(r$)
' Angle de vision (vertical ou horizontal) et distance de l'observateur
' permettent de définir l'emplacement de l'observateur
  repeat
    r$ = message_input$("Saisie de l'angle de vision vertical","Angle de vision vertical theta (-90° à 90°)"+nl$+"theta = ","40")
  until numeric(r$) = 1
  theta = val(r$)

  repeat
    r$ = message_input$("Saisie de l'angle de vision horizontal","Angle de vision horizontal phi (0° à 360°)"+nl$+"phi = ","350")
  until numeric(r$) = 1
  phi = val(r$)

  repeat
    r$ = message_input$("Saisie distance","Distance de l'observateur à l'origine ( > 5000 )"+nl$+"ray = ","10000")
  until numeric(r$) = 1
  ray = val(r$)

  ct = cos(theta*coef) : st = sin(theta*coef)
  cp = cos(phi*coef)  : sp = sin(phi*coef)
  xobs = ray*ct*cp    : yobs = ray*ct*sp  : zobs = ray*st

' Types du tracé :
' 1- tracé rapide: le paysage apparaitra en tracé fil-de-fer,sans élimination des parties cachées
' 2- parties cachées : même type de tracé, mais avec élimination des parties cachées
' 3- Surfaces éclairées : le programme demande la position de la source lumineuse
  repeat
    r$ = message_input$("Saisie type de tracé","Type de tracé (1, 2, 3)"+nl$+"1. Rapide"+nl$+"2. Parties cachées"+nl$+"3. Surface éclairée"+nl$+"type = ","3")
  until numeric(r$) = 1
  type = val(r$)

  if type = 3
    repeat
      r$ = message_input$("Saisie de l'angle vertical d'éclairage","Angle vertical d'éclairage alfa (-90° à 90°)"+nl$+"alfa = ","50")
    until numeric(r$) = 1
    alfa = val(r$)*coef
    repeat
      r$ = message_input$("Saisie de l'angle horizontal d'éclairage","Angle horizontal d'éclairage beta (0° à 360°)"+nl$+"beta = ","250")
    until numeric(r$) = 1
    beta = val(r$)*coef
    repeat
      r$ = message_input$("Saisie de la source lumineuse","Distance de la source lumineuse à l'origine ( > 5000)"+nl$+"r = ","10000")
    until numeric(r$) = 1
    r = val(r$)
    xecl = r*cos(alfa)*cos(beta)
    yecl = r*cos(alfa)*sin(beta)
    zecl = r*sin(alfa)
    repeat
      r$ = message_input$("Saisie des contours des facettes","Désirez-vous les contours des facettes (O/N)"+nl$+"Bord$ = ","O")
    until upper$(r$) = "O" or upper$(r$) = "N"
    bord$ = r$
  end_if
  caption 0,titre$ + "  **** !!!  VEUILLEZ    PATIENTER  !!!  **** " + "  <ESC>  POUR ARRETER"
return
' ******************************************************************************
generer_matrice:
  dn = power(2,mail)+1 : pas = dn - 1 : ech = 4000
  dim noeud(dn,dn)
  while pas > 1
    for i = 1 to dn-pas step pas
      for j = 1 to dn-i-pas+1 step pas
        a = i+pas/2 : b = j+pas/2 : cc = i+pas : d = j+pas : gosub calc_altitude
        noeud(i,b) = (noeud(i,j)+noeud(i,d))/2+alt        : gosub calc_altitude
        noeud(a,j) = (noeud(i,j)+noeud(cc,j))/2+alt        : gosub calc_altitude
        noeud(a,b) = (noeud(cc,j)+noeud(i,d))/2+alt
      next j
    next i
    pas = pas/2 : ech = ech/2
  end_while
  return
' ******************************************************************************
' *** Calcul aléatoire des altitudes
calc_altitude:
  alt = rnd(1)*ech
  if rnd(1) > profil/100 then alt = 0-alt
return
' ******************************************************************************
' *** Calcul du centrage
calc_centrage:

  xemin = 1000 : yemin = 1000 : xemax = 0-1000 : yemax = 0-1000
  if mail > 5 then stp = power(2,(mail-5)) : else : stp = 1
  for i = 1 to dn step stp
    for j = 1 to dn-i+1 step stp
      ci = i : cj = j
      gosub calc_coord_ecran
      if xxe < xemin then xemin = xxe
      if xxe > xemax then xemax = xxe
      if yye < yemin then yemin = yye
      if yye > yemax then yemax = yye
    next j
  next i

  ecrx = 990 : ecry = 690  : ' Dimension de l'écran graphique
  rap = 1 : ' 2.62 : ' permet de tenir compte de la déformation d'écran
  echx = ecrx /(xemax-xemin)/rap
  echy = ecry /(yemax-yemin)
  if echy < echx then echx = echy
  echy = echx : echx = echx*rap
  xcent = (ecrx-echx*(xemax+xemin))/2
  ycent = (ecry+10-echy*(yemax+yemin))/2
return
' ******************************************************************************
affichage:
  if type = 1 then gosub calc3 : return
  if phi >= 120 and phi < 240 then gosub test_phi :return
  if phi < 120 then gosub calc5 : return
  gosub calc4
return
' ******************************************************************************
' *** Calcul des coordonnées d'écran
calc_coord_ecran:
  gosub calc_coord_spatiales
coord_ecran:
  d = c(1)*cp*ct+c(2)*sp*ct+c(3)*st-ray
  xxe = (c(1)*sp-c(2)*cp)/d : yye = (c(1)*cp*st+c(2)*sp*st-c(3)*ct)/d
  yye = 0-yye
return
' ******************************************************************************
' *** Calcul des coordonnées spatiales
calc_coord_spatiales:
  c(1) = ((1-ci)/(dn-1)+1/3)*4000*sqr(3)
  c(2) = ((cj-1)+(ci-1)/2)*8000/(dn-1)-4000
  c(3) = noeud(ci,cj)
  if c(3) < mer then c(3) = mer
return
' ******************************************************************************
test_phi:
  if phi < 180 then deb = 1 : sens = 1 : else : fin = 1 : sens = 0-1
  for i = 2 to dn
    tst = 0
    if phi < 180 then fin = dn-i+1: else : deb = dn-i+1
    for j = deb to fin step sens
      if phi < 180 then gosub calc1: gosub calc2 : else : gosub calc2: gosub calc1
    next j
  next i
return
' ******************************************************************************
calc1:
  ci = i-1 : cj = j : gosub calc_coord_spatiales
  for l = 1 to 3
    v1(l) = c(l): pt(1,l) = c(l)
  next l
  ci = i : gosub calc_coord_spatiales
  for l = 1 to 3
    v1(l) = v1(l)-c(l): pt(2,l) = c(l)
  next l
  if tst = 0 then tst = 1 : return
  if phi < 180 then cj = j-1 : else : cj = j+1 : ci = i-1
  gosub calc_coord_spatiales
  for l = 1 to 3
    pt(3,l) = c(l)
  next l
  gosub remplis_facettes
return
' ******************************************************************************
calc2:
ci = i-1 : cj = j+1 : gosub calc_coord_spatiales
for l = 1 to 3
  v2(l) = c(l) : pt(1,l) = c(l)
next l
ci = i : cj = j : gosub calc_coord_spatiales
for l = 1 to 3
  v2(l) = v2(l)-c(l) : pt(2,l) = c(l)
next l
if tst = 0 then tst = 1 : return
if phi < 180 then ci = i-1 : else : cj =j+1
gosub calc_coord_spatiales
for l = 1 to 3
  pt(3,l) = c(l)
next l
gosub remplis_facettes
return
' ******************************************************************************
calc3:
  cc = 1
  for  i = 2 to dn
    for j = 1 to dn-i+1
      ci = i : cj = cj : gosub calc_coord_ecran
      xe(1) = xxe*echx+xcent : ye(1) = yye*echy+ycent
      ci = i-1 : gosub calc_coord_ecran
      xe(2) = xxe*echx+xcent : ye(2) = yye*echy+ycent
      cj = j+1 : gosub calc_coord_ecran
      xe(3) = xxe*echx+xcent : ye(3) = yye*echy+ycent
      gosub trace_contours
    next j
  next i
return
' ******************************************************************************
calc4:
  if phi < 300 then deb = 1 : sens = 1 : else : fin = 1 : sens = 0-1
  for i = dn-1 to 1 step -1
    tst = 0 : if phi < 300 then fin = i : else : deb = i
    for j = deb to fin step sens
      k = i+1-j
      if phi < 300 then gosub calc6 : gosub calc7 : else : gosub calc7 : gosub calc6
    next j
  next i
return
' ******************************************************************************
' ******************************************************************************
calc5:
  if phi > 60 then deb = 1 : sens = 1 : else : fin = 1 : sens = 0-1
  for i = 2 to dn
    tst = 0 : if phi > 60 then fin = dn-i+1 : else : deb = dn-i+1
    for j = deb to fin step sens
      if phi > 60 then gosub calc9 : gosub calc8 : else : gosub calc8 : gosub calc9
    next j
  next i
return
' ******************************************************************************
calc6:
  ci = j : cj = k + 1 : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = c(l) : pt(1,l) = c(l) : next l
  cj = k : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = v1(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi < 300 then ci = j-1 : cj = k+1 : else : ci = j+1
  gosub calc_coord_spatiales
  for l = 1 to 3 : pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
calc7:
  ci = j+1 : cj = k : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = c(l) : pt(1,l) = c(l) : next l
  ci = j : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = v2(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi < 300 then cj = k+1 : else : ci = j+1 : cj = k-1
  gosub calc_coord_spatiales
  for l = 1 to 3 : pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
calc8:
  ci = j+1 : cj = i-1 : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = c(l) : pt(1,l) = c(l) : next l
  ci = j : cj = i : gosub calc_coord_spatiales
  for l = 1 to 3 : v1(l) = v1(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi > 60 then cj = i-1 : else : ci = j+1
  gosub calc_coord_spatiales
  for l = 1 to 3 :  pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
calc9:
  ci = j : cj = i-1 : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = c(l) : pt(1,l) = c(l) : next l
  cj = i : gosub calc_coord_spatiales
  for l = 1 to 3 : v2(l) = v2(l)-c(l) : pt(2,l) = c(l) : next l
  if tst = 0 then tst = 1 : return
  if phi > 60 then ci = j-1 : else : ci = j+1 : cj = i-1
  gosub calc_coord_spatiales
  for l = 1 to 3 : pt(3,l) = c(l) : next l
  gosub remplis_facettes
return
' ******************************************************************************
' *** Tracé des contours
trace_contours:
  if cc = 1 then 2d_pen_color 0,0,0 : else : 2d_pen_color  xe(1) / 4,255,255 - (ye(1)/3)
  2d_line xe(1),ye(1),xe(2),ye(2)
  2d_line xe(2),ye(2),xe(3),ye(3)
  2d_line xe(3),ye(3),xe(1),ye(1)
return
' ******************************************************************************
' *** Calcul des couleurs et remplissagedes facettes
remplis_facettes:
  indm = 0 : ncoul = 0
  if pt(1,3) = mer and pt(2,3) = mer and pt(3,3) = mer then indm = 1 : goto suite1
  vn(1) = v1(2)*v2(3)-v1(3)*v2(2)
  vn(2) = v1(3)*v2(1)-v1(1)*v2(3)
  vn(3) = v1(1)*v2(2)-v1(2)*v2(1)
  vvn = sqr(vn(1)*vn(1)+vn(2)*vn(2)+vn(3)*vn(3))
  for n = 1 to 3 : bary(n) = (pt(1,n)+pt(2,n)+pt(3,n))/3 : next n
  vobs(1) = bary(1)-xobs : vobs(2) = bary(2)-yobs: vobs(3) = bary(3)-zobs
  prosc = vn(1)*vobs(1)+vn(2)*vobs(2)+vn(3)*vobs(3)
  vvobs = sqr(vobs(1)*vobs(1)+vobs(2)*vobs(2)+vobs(3)*vobs(3))
  cosang = prosc/(vvn*vvobs)
  ang = atn(sqr(1-cosang*cosang)/cosang)
  if ang > 0 then ncoul = 1
  if ncoul = 1 or type = 2  then goto suite1
  vecl(1) = bary(1)-xecl : vecl(2) = bary(2)-yecl : vecl(3) = bary(3)-zecl
  prosc = vn(1)*vecl(1)+vn(2)*vecl(2)+vn(3)*vecl(3)
  vvecl = sqr(vecl(1)*vecl(1)+vecl(2)*vecl(2)+vecl(3)*vecl(3))
  cosang = prosc/(vvn*vvecl)
  ang = atn(sqr(1-cosang*cosang)/cosang)
  if ang < 0 then ang = ang + pi
  coul = int(ang/fang+.5)
  if coul <> 0 then freq = 15/coul : else : freq = power(10,30)
suite1:
  for n = 1 to 3
    c(1) = pt(n,1) : c(2) = pt(n,2) : c(3) = pt(n,3)
    gosub coord_ecran
    xe(n) = int(xxe * echx + xcent + .5) : ye(n) = int(yye * echy + ycent + .5)
  next n
  if ye(2) >= ye(1) and ye(2) >= ye(3)
    exx = xe(1) : exy = ye(1) : xe(1) = xe(2) : ye(1) = ye(2)
    xe(2) = exx : ye(2) = exy : goto suite2
  end_if
  if ye(3) >= ye(1) and ye(3) >= ye(2)
    exx = xe(1) : exy = ye(1) : xe(1) = xe(3) : ye(1) = ye(3)
    xe(3) = exx : ye(3) = exy
  end_if
suite2:
  if ye(3) > ye(2)
    exx = xe(2) : exy = ye(2) : xe(2) = xe(3)
    ye(2) = ye(3) : xe(3) = exx : ye(3)= exy
  end_if
  ex1 = xe(1)-xe(2) : ey1 = ye(1)-ye(2)
  ex2 = xe(1)-xe(3) : ey2 = ye(1)-ye(3)
  ex3 = xe(2)-xe(3) : ey3 = ye(2)-ye(3)
  if ey1 <> 0 then stp1 = ex1/ey1
  if ey2 <> 0 then stp2 = ex2/ey2
  if ey3 <> 0 then stp3 = ex3/ey3
  pch = int(rnd(1)*freq+1.5) : comp = 1
  cpt2 = 0
  for cpt1 = 0 to ey1
    x1 = int(xe(1)-cpt1*stp1+.5) : x2 = int(xe(1)-cpt2*stp2+.5): y1 = ye(1)-cpt1
    if ey1 = 0 then x1 = xe(2)
    if ey2 = 0 then x2 = xe(3)
    gosub coloriage
  next cpt1
  for cpt1 = 1 to ey3
    x1 = int(xe(2)-cpt1*stp3+.5) : x2 = int(xe(1)-cpt2*stp2+.5): y1 = ye(1)-cpt2
    gosub coloriage
  next cpt1
  if (type=2 and ncoul=0 and indm=0)or(type=3 and ncoul = 1)or upper$(bord$)="O"
    cc = 1 : gosub trace_contours
  end_if
return

' ******************************************************************************
coloriage:
  if scancode = 27  then caption 0,"Arrêté par l'utilisateur" : end
  if (ncoul = 1 and type = 2) or (mod(y1,2) = 0 and indm = 1)
    cc = 1 : else : cc = 0
  end_if
  if type = 2 or indm = 1 or ncoul = 1
    if cc = 1 then 2d_pen_color 0,0,0 : else : 2d_pen_color xe(1) / 4,255,255 - (ye(1)/3)
    2d_pen_color xe(1) / 4,255,255 - (ye(1)/3)
    2d_line x1,y1,x2,y1
    cpt2 = cpt2 +1 : return
  end_if
  if x1 = x2
    if comp = int(pch+.5)
      2d_pen_color 0,0,0 : else : 2d_pen_color xe(1) / 4,255,255 - (ye(1)/3)
    end_if
    2d_point x1,y1
    cpt2 = cpt2 +1 : return
  end_if
for n = x1 to x2 step sgn(x2-x1)
    cc = 0
    if comp = int(pch+.5) then cc = 1 : pch = pch + freq
      if cc = 1 then 2d_pen_color 0,0,0 : else : 2d_pen_color xe(1) / 4,255,255 - (ye(1)/3)
    2d_point n,y1
    comp = comp + 1
  next n
  cpt2 = cpt2 + 1
return
' *************************************************************************
fin:
Caption 0,"Terminé"
end
' ******************************************************************************

...... 2d_pen_color mod(xe(1) / 4,255),255,mod(abs(255 - (ye(1)/3)),255)

 2d_pen_color mod(xe(1) / 4,255),rnd(255),mod(abs(255 - (ye(1)/3)),255)