#
# Wavefront & Shack-Hartmann Simulator
#
# Air bubbles & phase delays hypothesis
#
# Image (from original standalone program have been suppressed)
#  because plugin does not support binary extensions
#
# 1995/1996 by Laurent Demailly <dl@hplyot.obspm.fr>
#
# Thanks to Olivier Lai, Jean-Pierre Veran and Eric Gendron for physics hints
#

# note this was a quick hack primarly to create some
# figures for my thesis, the coding style is not not representative
# of my usual habits !

# Some TclX builtins replacements: (too slow...)
proc loop {var start end body} {
upvar $var v;
for {set v $start} {$v<$end} {incr v} {
uplevel $body;
}
}
# tcl-usage' faq random :
proc random {args} {
  global RNG_seed;
    
  set max 259200;
  set argcnt [llength $args];
  if { $argcnt < 1 || $argcnt > 2 } {
    error "wrong # args: random limit | seed ?seedval?"
  }
  if ![string compare [lindex $args 0] seed] {
    if { $argcnt == 2 } {
      set RNG_seed [lindex $args 1]
    } else {
      set RNG_seed [clock clicks] ; # poor...
    }
    return;
  }
  if ![info exists RNG_seed] {
    set RNG_seed [clock clicks] ; # poor...
  }
  set RNG_seed [expr ($RNG_seed*7141+54773) % $max]
  return [expr int(double($RNG_seed)*[lindex $args 0]/$max)]
}
proc lassign {list args} {
  set i 0;
  foreach vname $args {
    uplevel [list set $vname [lindex $list $i]];
    incr i;
  }
  lrange $list $i end
}

######### start of the program (french comments!) #########

# nombre de sous pupilles
set nb_pup 5
# nombre de point d'echantillons par sous pupilles
set nb_ptspup 3 ; #  9
# soit nombre de points :
#     +---+---+---+
#     123412341234X      <- nb_pup 4 , nb_ptspup 5
#         12345

# espace entre les points
set echx 20 ; # 10

# 'vitesse' dans milieu global (pas d'echantillonage en y)
set v0 20;
# vitesse lente dans bulle froide
set vlow 18;
# vitesse rapide dans bulle chaude
set vhigh 22;

# nombre de bulles
set nb_bulles 5 ; # 15

# sous echantillonage affichage: (affiche 1/$skipy calculs)
set skipy 1;

set offy 20;

# hauteur de la zone turbulente
set ty 220;

# generatoire aleatoire
random seed 45679;

# rayon minimal des bulles
set minr 10;
# ecart aleatoire maximum (maxr=minr+randr)
set randrstart 60 ;

# genere une nouvelle bulle 
proc bulle {} {
  global tx ty minr randr offy vlow vhigh;
# distribution 'kolmogorov' des tailles 
  set r [expr int($minr+pow([random [expr int(pow($randr,3))]],.333333))];
  set x [expr [random [expr $tx+2*($r-1)]]-$r+1];
  set y [expr $offy+$r+[random [expr $ty-(2*$r)]]];
  set v [expr [random 2]?$vlow:$vhigh];
  list $x $y $r $v;
}

# est-ce que la bulle (x,y,r) intercepte une bulle deja existante ?
proc isin {x y r} {
  global blist;
  foreach bul $blist {
    lassign $bul xb yb rb vb;
    set dx [expr $xb-$x];
    set dy [expr $yb-$y];
    set dist [expr $dx*$dx+$dy*$dy];
    if {$dist<($r+$rb)*($r+$rb)} {
       return 1
    }
  }
  return 0;
}

# vitesse du milieu au point (x,y)
proc getv {x y} {
global blist v0;
foreach bul $blist {
 lassign $bul xb yb rb vb;
 set dx [expr $xb-$x];
 set dy [expr $yb-$y];
 if {($dx*$dx+$dy*$dy)<=($rb*$rb)} {return $vb}
}
return $v0;
}

# nom du canvas de dessin
frame .tc
pack .tc -side bottom  -fill both -expand 1 -side right ;
set c ".tc.c"

#
# Construction de l'interface utilisateur
# permettant de changer les parametres, lancer et arretter la simu,...
#
frame .fl -relief groove -borderwidth 3;
label .fl.l1 -text "WF Simulation (toy)" ;
label .fl.l2 -text "Laurent Demailly"

button .pr -text "Start simulation" -command "do"

pack .fl .fl.l1 .fl.l2  -pady 1m -padx .5m;

# variables controlables et nom:

# with current tk_plugin, no toplevel so we have to fix
# dimensions...
# {nb_ptspup "Number of Samples per aperture"}
# {echx      "X Scale"}
# {nb_pup    "Number of Apertures"}
# {skipy     "Skip"}
# {offy      "Offy"}
# {ty        "Turbulence layer height"}

set lst {
 {nb_bulles "Number of Bubbles"}
 {v0        "Nominal Speed"}
 {vlow      "Speed in cold bubbles"}
 {vhigh     "Speed in hot (red) bubbles"}
 {minr      "Minimal bubble radius"}
 {randrstart "Aditional random radius"}
}

frame .fv;
pack .fv -padx .5m -pady .5m -fill x -expand 1;

# could be replaced by newer & nice grid manager... [this is old code of mine]
foreach el $lst {
 lassign $el nvar label;
 set f .fv.f$nvar;
 frame $f;
 label $f.l$nvar -text "$label :" -anchor w;
 entry $f.e$nvar -width 4 -textvariable $nvar -relief sunken ;
 pack $f.l$nvar -fill x -expand y -side left;
 pack $f.e$nvar -fill x -side left;
 pack $f -padx .5m -pady .5m -fill x -expand 1;
 }


pack .pr -pady .5m -padx .5m;

set font  "-*-helvetica-medium-r-normal--8-*"
set fontmap($font) [list "Helvetica-Bold" 8];
set colormap(red) "1 .3 .3 setrgbcolor"
set colormap(blue) "0 0 .4 setrgbcolor"

# etat de base (simu arrete'e)
proc reset {} {
.pr configure -state normal
catch {destroy .ps}
catch {destroy .lrem}
}

#
proc mk_canvas {} {
uplevel #0 {
set nb_pt [expr ($nb_ptspup*$nb_pup)+1]
set tx [expr $echx*($nb_pt-1)];

set toty [expr 3*$offy+$ty+150];
set totx [expr $tx+1];

catch {destroy $c}
canvas $c -width $totx -height $toty;
pack $c -fill y -expand 1;
$c create rectangle 0 0 $totx $toty -fill white -outline white;

# limites de la couche turbulente
$c create line 0 $offy $totx $offy -fill green -width 2;
$c create line 0 [expr $offy+$ty] $totx [expr $offy+$ty] -fill green -width 2;
}
}

mk_canvas;

#
# Proc do: lance effectivement la simu
#
proc do {} {
uplevel #0 {
set stop 0;
.pr configure -state disabled

set status "Init...";
label .lrem -textvariable status;
pack .lrem -padx .5m -pady .5m;

button .ps -text Stop -command "set stop 1"
pack .ps -padx .5m -pady .5m;
update;
mk_canvas;
#
#
#
#
# Generation des bulles
#
set blist {};
set randr $randrstart;
loop i 0 $nb_bulles {
  set status "Remaining [expr $nb_bulles-$i] bubble to place...";
  set j 0;
  while 1 {
    set bul [bulle];
    lassign $bul x y r v;
    set b [$c create oval [expr $x-$r] [expr $y-$r] [expr $x+$r] [expr $y+$r] -width 3 -outline pink];
    if {$stop} {reset;return;}
    if {![isin $x $y $r]} break;
    update
    $c delete $b;
    incr j;
    if {($j>3) && ($randr>5)} {incr randr -1}
  }
  if {$v>$v0} {set col red} else {set col blue};
  $c itemconfig $b -outline $col;
  lappend blist $bul;
  update;
}

set status "Generating plane wavefront..."
update;

#puts "working on $nb_pt points";
set ptlist {};
# starting point
loop i 0 $nb_pt {
  lappend ptlist [list [expr $echx*$i] $v0];
}

eval $c create line [join $ptlist " "];
# 
set i 0;

# Y du plan des pupilles
set yp [expr 3*$offy+$ty];


set status "Propagating WF ";
#
# de'placement du front d'onde :
#
while 1 {
  incr i;
  set newlist {};
  set ymin 9999;
  foreach pt $ptlist {
   lassign $pt x y;
   set ny [expr $y+[getv $x $y]];
   if {$ny<$ymin} {set ymin $ny}
   lappend newlist "$x $ny";
  }
  if {$stop} {reset;return;}
  set ptlist $newlist;
  if {$i%$skipy==0} {
    eval $c create line [join $ptlist " "];
    append status ".";
    update;
    if {$ymin>=$yp} break;
  }
}

#puts "done $i iters";
set status "WF analysis & correction...";
update;

# dessine les micro lentilles :
set szpup [expr $echx*$nb_ptspup];
set szpup2 [expr $szpup/2];
set cx 0;
set dy 6;
set yp1 [expr $yp-$dy]
set yp2 [expr $yp+$dy]
$c create rectangle 0 $yp $totx $toty -fill white -outline white;

lassign [lindex $ptlist 0] x0 y0;
set mlist {}

loop i 0 $nb_pup {
  $c create oval $cx $yp1 [expr $cx+$szpup] $yp2 -width 2 -outline brown -fill white;
  $c create line $cx [expr $yp+100] $cx [expr $yp-$offy/2] -fill brown;
  set cx2 [expr $cx+$szpup2];
  set p1 [lindex $ptlist [expr $i*$nb_ptspup]];
  set p2 [lindex $ptlist [expr ($i+1)*$nb_ptspup]];
#  puts "$p1 $p2";
  lassign $p1 x1 y1;
  lassign $p2 x2 y2;
  # calcul du miroir segmente'
  lappend mlist "$x1 [expr ($y1-$y0)/2.]"
  lappend mlist "$x2 [expr ($y2-$y0)/2.]"
  # vecteur normal ( $y1-$y2 , $x2-$x1 )
  $c create line $cx2 $yp [expr $cx2+.5*($y1-$y2)] [expr $yp+.5*($x2-$x1)];
  incr cx $szpup;
}
  $c create line $cx [expr $yp+100] $cx [expr $yp-$offy/2] -fill brown;
  update
  lassign [correct $nb_pup $nb_ptspup $ptlist] nptlist ptcorr;
  set n [eval $c create line [join $mlist " "] -fill black -width 2];
  $c move $n 0 [expr $yp+140];
  set n [eval $c create polygon [join $mlist " "] $cx $offy 0 $offy -fill orange1];
  $c move $n 0 [expr $yp+140];
  set n [eval $c create line [join $nptlist " "] -fill darkgreen];
  $c move $n 0 [expr $yp+100];
  set n [eval $c create line [join $ptcorr " "] -fill gold];
  $c move $n 0 [expr $yp+100];
  catch {$c create text [expr $totx-20] [expr $toty-5] -font $font -text "(c)DL"}
}
reset
}

# en cas d'erreur...
#proc tkerror {info} {
#reset;
#tkerror.tk $info;
#}


set lscale 20.;
set lambda 0.5;


#
# "corrige" le front d'onde : mirroir segmenté.
# et retire un pseudo piston (moyenne des extremitées) de pt source
proc correct {nb_pup nb_ptspup ptlist} {
 set ptlr {};
 set ptloc {};
 set n0 0;
 set n1 $nb_ptspup;
 lassign [lindex $ptlist 0] x0 y0;
 lassign [lindex $ptlist [expr $nb_pup*$nb_ptspup]] xn yn;
 set cor [expr ($y0+$yn)/2];
 loop i 0 $nb_pup {
   lassign [lindex $ptlist $n0] x1 y1;
   lassign [lindex $ptlist $n1] x2 y2;
   lappend ptlr "$x1 0."
   lappend ptloc "$x1 [expr $y1-$cor]";
   set pente [expr 1.*($y2-$y1)/($x2-$x1)]
   loop j 1 $nb_ptspup {
     lassign [lindex $ptlist [expr $n0+$j]] x y;
     lappend ptlr "$x [expr $y-$y1-($x-$x1)*$pente]"
     lappend ptloc "$x [expr $y-$cor]";
   }
   incr n0 $nb_ptspup;
   incr n1 $nb_ptspup;
 }
 lappend ptlr "$x2 0.";
 lappend ptloc "$x2 [expr $y2-$cor]";
 list $ptloc $ptlr;
}