/* parametres */
dist_gaz_gaz := 1.0 ;; /* relation de voisinage */
dist_mem_mem := 1.0 ;;
dist_mem_gaz := 2.5 ;;

param_electrom := 5.0 ;; /* puissance des forces */ 
param_chimique := (-50.0) ;;

nb_gaz := 30 ;;
angle_membranes := 0.12 ;;

nb_pas := 2000 ;;

/* Déclaration des types */

/*record Particule = { px:float , py:float , vx:float , vy:float , m:float } ;;
record Membrane  = Particule + { membrane , ~gaz , idx:int , v_idx1:int , v_idx2:int } ;;
record Gaz       = Particule + { ~membrane , gaz } ;;
record Forces    = { fx:float , fy:float } ;;*/

record Particule = { px , py , vx , vy , m } ;;
record Membrane  = Particule + { membrane , ~gaz , idx , v_idx1 , v_idx2 } ;;
record Gaz       = Particule + { ~membrane , gaz } ;;
record Forces    = { fx , fy } ;;

delaunay(2,true) E2 =
  \e.if Particule(e)
     then (e.px, e.py)
     else ?("bad element type for E2 delaunay type")
     fi
;;

fun dist (p1,p2) = sqrt( (p1.px - p2.px)*(p1.px - p2.px) + (p1.py - p2.py)*(p1.py - p2.py) ) ;;
fun dist2(p1,p2) = (p1.px - p2.px)*(p1.px - p2.px) + (p1.py - p2.py)*(p1.py - p2.py) ;;
fun dist4(p1,p2) = let d = dist2(p1,p2) in d*d ;;

fun prox(p1,p2) = 
 (
  if (Membrane(p1) & Membrane(p2)) then
    ((p2.idx == p1.v_idx1) | (p2.idx == p1.v_idx2)| (dist(p1,p2) < dist_mem_mem) )
  else
    if (Gaz(p1) & Gaz(p2)) then
      (dist(p1,p2) < dist_gaz_gaz)
    else
      (dist(p1,p2) < dist_mem_gaz)
    fi
  fi
);;

proximal P2 = prox ;;

/* Initialisation d'une collection test */

fun new_gaz_part() = (
  { px = random(10.0) - 5.0 , py = random(10.0) - 5.0 , vx = 0.0 , vy = 0.0 , m = 2.0 , gaz }
);;

fun cloud_gaz(n,acc) = let rayon = 3.5 in
    (
  if (n==0) then
    acc
  else
    let np = new_gaz_part() in
      if ((sqrt(np.px*np.px + np.py*np.py)) < rayon) then
	np , cloud_gaz(n-1,acc)
      else
	cloud_gaz(n,acc)
      fi
  fi
);;

fun new_membrane_parts(r,t,i) = (
  let eps = 0.0 in//random(0.002) - 0.001 in
  let c = r * cos(t) + eps in
  let s = r * sin(t) - eps in
    { px =  c  , py =  s  , vx = 0.0 , vy = 0.0 , m = 10.0 , membrane , idx = i }
);;

fun make_membrane(r,t,i) = (
  if (t >= (2*PI)) then
    seq:()
  else (
    nb_part := nb_part + 1 ;
    new_membrane_parts(r,t,i) , make_membrane(r,t+angle_membranes,i+1)
  ) fi
);;

trans init_membrane [n] = {
  p => p + { v_idx1 = (let x=((p.idx-1) % n) in if (x == -1) then (n-1) else x fi), v_idx2 = ((p.idx+1) % n) }
};;

trans incr [n] = {
  p => p+ {idx = p.idx +n ,v_idx1 = p.v_idx1 +n , v_idx2 = p.v_idx2 +n }
} ;;
trans translation [x,y] = {
  p => p+ {px = p.px +x ,py = p.py +y }
} ;;
trans vitesse [x,y] = {
  p => p+ {vx = p.vx +x ,vy = p.vy +y }
} ;;

fun build_mem(s) = (
  if (s == seq:()) then
    P2:()
  else
    s.(0),(build_mem(tl(s)))
  fi
);;


nb_part := 0 ;;
init := make_membrane(4.0,0.0,0) ;;
init := init_membrane[n=nb_part](init) ;;
init := cloud_gaz(nb_gaz,init);;
init := vitesse[x=1.0,y=0.0](init) ;;

nb_part := 0 ;;
init2 := make_membrane(4.0,0,0) ;;
init2 := init_membrane[n=nb_part](init2) ;;
init2 := incr[n=200](init2) ;;
init2 := cloud_gaz(nb_gaz,init2);;
init2 := translation[x=9,y=4](init2) ;;
init2 := vitesse[x=-1.0,y=-0.0](init2) ;;

init := init,init2 ;;

stdout << "debut de proxymification\n";;
init := build_mem(init);;
stdout << "fin de proxymification\n";;

//init2 := cloud_gaz(nb_gaz,build_mem(init2));;

//init := init2,init ;;

/*init := 
{ gaz = <undef>, m = 2.000000, px = 7.0, py = 0.0, vx = -1.5, vy = 0.0},
{ idx = 103, m = 10.000000, membrane = <undef>, px = 8.0, py = 0.0, v_idx1 = 102, v_idx2 = 100, vx = -1.5, vy = 0.0},
{ idx = 102, m = 10.000000, membrane = <undef>, px = 7.0, py = 1.0, v_idx1 = 101, v_idx2 = 100, vx = -1.5, vy = 0.0},
{ idx = 101, m = 10.000000, membrane = <undef>, px = 6.0, py = 0.0, v_idx1 = 100, v_idx2 = 102, vx = -1.5, vy = 0.0},
{ idx = 100, m = 10.000000, membrane = <undef>, px = 7.0, py = -1.0, v_idx1 = 103, v_idx2 = 101, vx = -1.5, vy = 0.0},init;;
*/

/*
init :=
{ gaz = <undef>, m = 2.000000, px = -0.347404, py = 0.244785, vx = 0.077590, vy = 0.100764 },
{ gaz = <undef>, m = 2.000000, px = 0.172376, py = -0.395843, vx = 0.073952, vy = 0.068311 },
{ gaz = <undef>, m = 2.000000, px = -0.000751, py = -0.332799, vx = -0.045003, vy = -0.065263 },
{ gaz = <undef>, m = 2.000000, px = -0.067660, py = -0.195102, vx = -0.076084, vy = -0.114854 },
{ gaz = <undef>, m = 2.000000, px = 0.057218, py = 0.165812, vx = 0.062945, vy = -0.168024 },
{ gaz = <undef>, m = 2.000000, px = 0.437530, py = -0.119466, vx = -0.117636, vy = -0.140812 },
{ gaz = <undef>, m = 2.000000, px = 0.107131, py = -0.135808, vx = -0.160444, vy = -0.173232 },
{ gaz = <undef>, m = 2.000000, px = -0.138949, py = 0.066157, vx = 0.096939, vy = 0.307549 },
{ gaz = <undef>, m = 2.000000, px = -0.346356, py = -0.041974, vx = -0.206451, vy = -0.104812 },
{ gaz = <undef>, m = 2.000000, px = 0.434552, py = 0.142285, vx = -0.145369, vy = -0.249434 },
{ gaz = <undef>, m = 2.000000, px = -0.160002, py = -0.126539, vx = 0.003523, vy = -0.223521 },
{ gaz = <undef>, m = 2.000000, px = 0.234164, py = -0.290308, vx = 0.137238, vy = 0.087824 },
{ gaz = <undef>, m = 2.000000, px = -0.253808, py = -0.204956, vx = -0.147315, vy = -0.202217 },
{ gaz = <undef>, m = 2.000000, px = -0.144243, py = 0.352594, vx = -0.117148, vy = -0.012795 },
{ gaz = <undef>, m = 2.000000, px = 0.046905, py = -0.001720, vx = 0.065489, vy = 0.065626 },
{ idx = 0, m = 10.000000, membrane = <undef>, px = 2.347031, py = 0.239015, v_idx1 = 25, v_idx2 = 1, vx = -0.067727, vy = 0.109896+0.2 },
{ idx = 1, m = 10.000000, membrane = <undef>, px = 2.367808, py = 0.716134, v_idx1 = 0, v_idx2 = 2, vx = -0.106739, vy = 0.002791+0.2 },
{ idx = 2, m = 10.000000, membrane = <undef>, px = 2.000285, py = 1.356133, v_idx1 = 1, v_idx2 = 3, vx = -0.003567, vy = 0.026580+0.2 },
{ idx = 3, m = 10.000000, membrane = <undef>, px = 1.548756, py = 1.785272, v_idx1 = 2, v_idx2 = 4, vx = -0.029194, vy = -0.039360 +0.2},
{ idx = 4, m = 10.000000, membrane = <undef>, px = 0.821233, py = 2.127665, v_idx1 = 3, v_idx2 = 5, vx = 0.067145, vy = -0.052441+0.2 },
{ idx = 5, m = 10.000000, membrane = <undef>, px = 0.246632, py = 2.364512, v_idx1 = 4, v_idx2 = 6, vx = -0.054890, vy = -0.028477+0.2 },
{ idx = 6, m = 10.000000, membrane = <undef>, px = -0.543679, py = 2.199272, v_idx1 = 5, v_idx2 = 7, vx = 0.126439, vy = 0.097490+0.2 },
{ idx = 7, m = 10.000000, membrane = <undef>, px = -1.070309, py = 2.084454, v_idx1 = 6, v_idx2 = 8, vx = -0.077596, vy = 0.032797+0.2 },
{ idx = 8, m = 10.000000, membrane = <undef>, px = -1.823991, py = 1.793278, v_idx1 = 7, v_idx2 = 9, vx = 0.159582, vy = -0.133411+0.2 },
{ idx = 9, m = 10.000000, membrane = <undef>, px = -2.025106, py = 1.349159, v_idx1 = 8, v_idx2 = 10, vx = -0.048207, vy = 0.185766+0.2 },
{ idx = 10, m = 10.000000, membrane = <undef>, px = -2.351005, py = 0.849306, v_idx1 = 9, v_idx2 = 11, vx = 0.000254, vy = 0.032464 +0.2},
{ idx = 11, m = 10.000000, membrane = <undef>, px = -2.332896, py = 0.379105, v_idx1 = 10, v_idx2 = 12, vx = -0.015776, vy = 0.071074+0.2 },
{ idx = 12, m = 10.000000, membrane = <undef>, px = -2.279949, py = -0.100430, v_idx1 = 11, v_idx2 = 13, vx = 0.078300, vy = 0.097275 +0.2},
{ idx = 13, m = 10.000000, membrane = <undef>, px = -2.294624, py = -0.559415, v_idx1 = 12, v_idx2 = 14, vx = -0.080657, vy = -0.015653 +0.2},
{ idx = 14, m = 10.000000, membrane = <undef>, px = -2.098781, py = -1.010770, v_idx1 = 13, v_idx2 = 15, vx = 0.059342, vy = -0.057670 +0.2},
{ idx = 15, m = 10.000000, membrane = <undef>, px = -1.899193, py = -1.508682, v_idx1 = 14, v_idx2 = 16, vx = -0.075261, vy = -0.066641 +0.2},
{ idx = 16, m = 10.000000, membrane = <undef>, px = -1.398293, py = -1.969579, v_idx1 = 15, v_idx2 = 17, vx = -0.068734, vy = -0.020676+0.2 },
{ idx = 17, m = 10.000000, membrane = <undef>, px = -1.158806, py = -2.207898, v_idx1 = 16, v_idx2 = 18, vx = 0.108011, vy = -0.024207+0.2 },
{ idx = 18, m = 10.000000, membrane = <undef>, px = -0.600661, py = -2.331795, v_idx1 = 17, v_idx2 = 19, vx = -0.149615, vy = 0.007197+0.2 },
{ idx = 19, m = 10.000000, membrane = <undef>, px = -0.274231, py = -2.292526, v_idx1 = 18, v_idx2 = 20, vx = 0.032829, vy = -0.256957+0.2 },
{ idx = 20, m = 10.000000, membrane = <undef>, px = 0.274390, py = -2.488804, v_idx1 = 19, v_idx2 = 21, vx = 0.013450, vy = 0.040196+0.2 },
{ idx = 21, m = 10.000000, membrane = <undef>, px = 0.803741, py = -2.379167, v_idx1 = 20, v_idx2 = 22, vx = 0.106300, vy = 0.017300+0.2 },
{ idx = 22, m = 10.000000, membrane = <undef>, px = 1.304850, py = -1.996203, v_idx1 = 21, v_idx2 = 23, vx = 0.009812, vy = 0.085522+0.2 },
{ idx = 23, m = 10.000000, membrane = <undef>, px = 1.655333, py = -1.701119, v_idx1 = 22, v_idx2 = 24, vx = 0.129464, vy = -0.004711+0.2 },
{ idx = 24, m = 10.000000, membrane = <undef>, px = 2.089812, py = -1.121035, v_idx1 = 23, v_idx2 = 25, vx = -0.066847, vy = -0.036875 +0.2},
{ idx = 25, m = 10.000000, membrane = <undef>, px = 2.374002, py = -0.474420, v_idx1 = 24, v_idx2 = 0, vx =0.006441, vy = 0.059752+0.2 },
{ idx = 100, m = 10.000000, membrane = <undef>, px = 4.0, py = 0.0, v_idx1 = 102, v_idx2 = 101, vx = -1.0, vy = -0.036875 +0.2},
{ idx = 101, m = 10.000000, membrane = <undef>, px = 5.0, py = 0.5, v_idx1 = 100, v_idx2 = 102, vx = -1.0, vy = -0.036875 +0.2},
{ idx = 102, m = 10.000000, membrane = <undef>, px = 6.0, py = -0.5, v_idx1 = 101, v_idx2 = 100, vx = -1.0, vy = -0.036875 +0.2},

 P2:() ;;
*/

trans init_speed = {
  p => p+{vx = p.vx+1.0}
};;

trans average_gaz [avr = ({ m = 0.0, px = 0.0, py = 0.0, vx = 0.0, vy = 0.0, n=0,
			    px_min = max_int, px_max = min_int,
			    py_min = max_int, py_max = min_int
			  }) ] = {
  g:Gaz / ( avr :=  
	    {
	      m  = avr.m  + g.m ,
	      px = avr.px + g.px ,
	      py = avr.py + g.py ,
	      vx = avr.vx + g.vx ,
	      vy = avr.vy + g.vy ,
	      n  = avr.n  + 1,
	      px_min = if (g.px < avr.px_min) then g.px else avr.px_min fi,
	      py_min = if (g.py < avr.py_min) then g.py else avr.py_min fi,
	      px_max = if (g.px > avr.px_max) then g.px else avr.px_max fi,
	      py_max = if (g.py > avr.py_max) then g.py else avr.py_max fi
	    } ; false ) => !!(false) ;
  _ =>
     let rayon = sqrt ( (avr.px_max - avr.px_min)*(avr.px_max - avr.px_min) +
			(avr.py_max - avr.py_min)*(avr.py_max - avr.py_min)) / 2 in
       Return(avr , {
		gaz ,
		m  = avr.m  ,
		px = avr.px / avr.n ,
		py = avr.py / avr.n ,
		vx = avr.vx / avr.n ,
		vy = avr.vy / avr.n
	      } , rayon )
};;

//average_gaz(init) ;;

//init := init_speed(init);;


fun addforces (f1,f2) =
     { fx = f1.fx + f2.fx ,
       fy = f1.fy + f2.fy } ;;  

/* ---------------------------------- */
/* Forces agissant sur les particules */
/* ---------------------------------- */


fun force_electrom(a,b) = 
    (let d2 = dist2(a,b) in
     { fx = param_electrom * (a.px - b.px) / (d2),
       fy = param_electrom * (a.py - b.py) / (d2) }) ;; 



fun force_chimique(a,b) = 
   
(let d = dist(a,b) in
//     { fx = param_chimique * (a.px - b.px) * (1- 1/ d),
//       fy = param_chimique * (a.py - b.py) * (1- 1/ d) }) ;; 
     { fx = param_chimique * (a.px - b.px) * d,
       fy = param_chimique * (a.py - b.py) * d }) ;; 


trans forces[t = 0.02] = {
      p:Gaz => (
	    let f = 
	    neighborsfold( 
	      (\p'.\r.(let f_elem = force_electrom(p,p') in addforces(f_elem,r))),
              {fx = 0.0, fy = 0.0}, 
              p )
	    in let v = {vx = p.vx+f.fx*t/p.m , vy = p.vy+f.fy*t/p.m } 
	    in p + { px = p.px+v.vx*t , py = p.py+v.vy*t } + v 
 );
       p:Membrane => (
	    let f = 
	    neighborsfold( 
	      (\p'.\r.(let f_elem = if Membrane(p')
                                    then (let voisid =  (p'.idx) in if ((voisid  == (p.v_idx1))  | (voisid == (p.v_idx2)) )
                                         then addforces (force_electrom(p,p'),force_chimique(p,p'))
                                         else force_electrom(p,p') fi) 
                                    else force_electrom(p,p') fi
                        in addforces(f_elem,r))),			     
              {fx = 0.0, fy = 0.0}, 
              p )
	    in let v = {vx = p.vx+f.fx*t/p.m , vy = p.vy+f.fy*t/p.m } 
	    in p + { px = p.px+v.vx*t , py = p.py+v.vy*t } + v 
 )	

};;

/*
                                         */
"OK";;

/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/


/* Fonction d'affichage d'un delaunay... */

fun show_node(p) = (
  "Translated { Translation <"
  + p.px + ", " + p.py + ", " + (0.0)
  + "> Geometry Sphere { Radius "
  + if (Membrane(p))
    then (0.1)
    else (0.1)
    fi
  + " Slices 5 Stacks 5 "
  + if (Membrane(p))
    then "Color <" + (1.0) + ", "+ (0.0) + ", " + (0.0) + "> } }"
    else "Color <" + (0.0) + ", "+ (1.0) + ", " + (0.0) + "> } }"
    fi
);;

trans draw_edge[acc] =
{
  x, y /
    (
      acc :=
       "\tPolyline { PointList [ <"+x.px+"  , "+x.py+"  , 0 > , <"+y.px+"  , "+y.py+"  , 0 > ] Color <" + (0.0) + ", "+ (0.0) + ", " + (0.0) + "> }\n"
       + acc ;
      false
    ) => (stdout << "ERROR : draw_edges\n" ; !!(false)) ;
  _ => Return(acc)
};;

fun show(f, c, i) = (
  uprint_coll(f, c, show_node,
              "Scene a"+i+" {\n\t",
              "\n\t",
              "\n"+
              (draw_edge[acc=""](c))+"\n"+
              "}\n\n")
) ;;

/* Sortie graphique avec Imoview */

out_file  := "./blob.imo" ;;
scene_num := 0 ;;

fun evol(c) = (
  show(out_file, c, scene_num);
  scene_num:= scene_num+1 ;
  stdout << scene_num;
  forces(c)
) ;;

fun fff(sn) = (
  if (sn==0) then
    out_file << "Replace { Show a0 }\n"
  else (
    fff(sn - 1);
    out_file << "Replace { Show a"+sn+" }\n"
  ) fi
);;

fun save_coll(f,c) = (
  uprint_coll(f,c,(\e.(e)),"init := \n",",\n",", P2:() ;;\n")
);;

res := evol['iter = nb_pas](init); 0;;
show(out_file, res, scene_num);
stdout << scene_num;
fff(scene_num);;

//save_coll("/tmp/init.mgs",res);;

//system("~/BIN/imoview.opt "+out_file);;