text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* One should be able to simplify this by evaluating explicitly the Stoll
   decomposition of the exchange functional *)

$define lda_c_pw_params
$define lda_c_pw_modified_params
$include "lda_c_pw.mpl"

b97mv_ux    := (mgamma, x) -> mgamma*x^2/(1 + mgamma*x^2):

(* article uses t = 2 tau convention *)
b97mv_wx_ss := (t, dummy) -> (K_FACTOR_C - t)/(K_FACTOR_C + t):
b97mv_wx_os := (ts0, ts1) -> (K_FACTOR_C*(ts0 + ts1) - 2*ts0*ts1)/(K_FACTOR_C*(ts0 + ts1) + 2*ts0*ts1):

b97mv_g := (mgamma, wx, cc, n, xs, ts0, ts1) ->
  add(cc[i][1]*wx(ts0, ts1)^cc[i][2]*b97mv_ux(mgamma, xs)^cc[i][3], i=1..n):

b97mv_fpar  := (rs, z, xs0, xs1, ts0, ts1) ->
  + lda_stoll_par(f_pw    , rs,  z,  1) * b97mv_g(b97mv_gamma_ss, b97mv_wx_ss, b97mv_par_ss, b97mv_par_n, xs0, ts0, 0)
  + lda_stoll_par(f_pw    , rs, -z, -1) * b97mv_g(b97mv_gamma_ss, b97mv_wx_ss, b97mv_par_ss, b97mv_par_n, xs1, ts1, 0):

b97mv_fos := (rs, z, xs0, xs1, ts0, ts1) ->
  lda_stoll_perp(f_pw, rs, z)
  * b97mv_g(b97mv_gamma_os, b97mv_wx_os, b97mv_par_os, b97mv_par_n, sqrt(xs0^2 + xs1^2)/sqrt(2), ts0, ts1):

b97mv_f :=  (rs, z, xs0, xs1, ts0, ts1) ->
  + b97mv_fpar(rs, z, xs0, xs1, ts0, ts1)
  + b97mv_fos(rs, z, xs0, xs1, ts0, ts1):


text:
`is_element/YC` := (N::posint) -> (A::set) -> proc(x)
 local a,b,i,j,k,t,G,J,E,B;
 global reason;

 if not(type(x,table)) then
  reason := [convert(procname,string),"x is not a table",x];
  return false;
 fi;

 if map(op,{indices(x)}) <> A then
  reason := [convert(procname,string),"x is not indexed by A",x,A];
  return false;
 fi;

 for a in A do
  if not `is_element/R`(N)(x[a]) then
   reason := [convert(procname,string),"x[a] is not in R^N",a,x[a],N];
   return false;
  fi;
 od;

 for i from 1 to N do
  if min(seq(x[a][i],a in A)) <> 0 then
   reason := [convert(procname,string),"min(x[i]) <> 0",i,x];
   return false;
  fi;
 od;

 G := `gaps/F`(N)(A)(x);

 if max(G) <> 1 then
  reason := [convert(procname,string),"max(G) <> 1",G];
  return false;
 fi;

 J := select(j -> G[j] = 1,[seq(k,k=1..N)]);
 E := {seq(seq([a,b],b in A),a in A)};
 for j in J do
  t := table();
  for a in A do t[a] := x[a][j]; od;
  B := blur_preorder(A)(t,1);
  E := E intersect B intersect `op/preord`(A)(B);
 od:

 if E <> {seq([a,a],a in A)} then
  reason := [convert(procname,string),"E <> Delta",E];
  return false;
 fi;

 return true;
end;

`is_equal/YC` := (N::posint) -> (A::set) -> proc(x,y)
 local a;
 global reason;

 for a in A do
  if x[a] <> y[a] then
   reason := [convert(procname,string),"x[a] <> y[a]",a,x[a],y[a]];
   return false;
  fi;
 od;

 return true;
end;

`inc/YC/F` := (N::posint) -> (A::set) -> (x) -> x;

`mu_bar/YC/PP` := (N::posint) -> (A::set) -> proc(x)
 local i,j,k,p,m,a,b,L,A0,T0,g0,g,t,s,AA,R,RR;

# DEBUG();

 L := NULL;
 m := nops(A);

 for i from 1 to N do
  A0 := sort([op(A)],(a,b) -> x[a][i] <= x[b][i]);
  T0 := [seq(x[A0[j]][i],j=1..m)];
  for j from 1 to m-1 do
   L := L,[i,A0[j],A0[j+1],T0[j+1]-T0[j]];
  od;
 od;

 L := sort([L],(a,b) -> a[4] <= b[4]);
 k := nops(L);

 R := table():
 t := table():

 R[1] := `mu/F/PP`(N)(A)(x);
 t[1] := L[1][4];
 s := table():
 s[R[1]] := t[1];

 AA := {seq(seq([a,b],b in A),a in A)};

 for p from 2 to k do
  i,a,b,g := op(L[p-1]);
  RR := R[p-1][i];
  RR := RR union select(xy -> member([xy[1],b],RR) and member([a,xy[2]],RR),AA);
  R[p] := [seq(R[p-1][j],j=1..i-1),RR,seq(R[p-1][j],j=i+1..N)];
  t[p] := L[p][4] - L[p-1][4];
  if t[p] > 0 then s[R[p]] := t[p]; fi;
 od;

 return eval(s);
end:

`sigma_bar/PP/YC` := (N::posint) -> (A::set) -> proc(t)
 local x,y,a,C,s,R;

 x := table();
 for a in A do x[a] := [0$N]; od;

 C := map(op,[indices(t)]);
 for R in C do
  s := t[R];
  y := `sigma/PP/F`(N)(A)(R);
  for a in A do
   x[a] := x[a] +~ s *~ y[a];
  od;
 od;

 return eval(x);
end:

ext: mpl


text:
TITLE
    Invest4Scen;

INDEX
    invtype := (Stock, Bond) -> (S,B);

STOCHASTIC

STAGES
    period := 1..4;

SCENARIOS
    sc := 1..8;

TREE
    BINARY;

PROBABILITIES
    prob[sc] := ALLEQUAL;

RANDOM DATA
    Returns[invtype, period, sc] :=
      [Stock, 1, 1, 1.25,
       Stock, 1, 2, 1.25,
       Stock, 1, 3, 1.25,
       Stock, 1, 4, 1.25,
       Stock, 1, 5, 1.06,
       Stock, 1, 6, 1.06,
       Stock, 1, 7, 1.06,
       Stock, 1, 8, 1.06,
       Stock, 2, 1, 1.25,
       Stock, 2, 2, 1.25,
       Stock, 2, 3, 1.06,
       Stock, 2, 4, 1.06,
       Stock, 2, 5, 1.25,
       Stock, 2, 6, 1.25,
       Stock, 2, 7, 1.06,
       Stock, 2, 8, 1.06,
       Stock, 3, 1, 1.25,
       Stock, 3, 2, 1.06,
       Stock, 3, 3, 1.25,
       Stock, 3, 4, 1.06,
       Stock, 3, 5, 1.25,
       Stock, 3, 6, 1.06,
       Stock, 3, 7, 1.25,
       Stock, 3, 8, 1.06,
       Bond, 1, 1, 1.14,
       Bond, 1, 2, 1.14,
       Bond, 1, 3, 1.14,
       Bond, 1, 4, 1.14,
       Bond, 1, 5, 1.12,
       Bond, 1, 6, 1.12,
       Bond, 1, 7, 1.12,
       Bond, 1, 8, 1.12,
       Bond, 2, 1, 1.14,
       Bond, 2, 2, 1.14,
       Bond, 2, 3, 1.12,
       Bond, 2, 4, 1.12,
       Bond, 2, 5, 1.14,
       Bond, 2, 6, 1.14,
       Bond, 2, 7, 1.12,
       Bond, 2, 8, 1.12,
       Bond, 3, 1, 1.14,
       Bond, 3, 2, 1.12,
       Bond, 3, 3, 1.14,
       Bond, 3, 4, 1.12,
       Bond, 3, 5, 1.14,
       Bond, 3, 6, 1.12,
       Bond, 3, 7, 1.14,
       Bond, 3, 8, 1.12];

DATA
    InitialWealth = 55;
    Goal := 80;
    SurplusReward := 1;
    ShortPenalty  := 4;

VARIABLES
    Invest[period<LAST(period),invtype] -> x;

LASTSTAGE VARIABLES
    Shortage -> Sh;
    Excess -> Ex;

MODEL

   MAX obj = SurplusReward*Excess - ShortPenalty*Shortage;

SUBJECT TO

   InitInvest[period=1] -> InitInv:
       InitialWealth
     =
       SUM(invtype: Invest);

   InvestBal[period=2..LAST(period)-1] -> InvBal:
       SUM(invtype: Invest[period-1] * Returns[period-1])
     =
       SUM(invtype: Invest);

   EndInvest[period=LAST(period)] -> EndInv:
       SUM(invtype: Invest[period-1]* Returns[period-1])  + Shortage - Excess
     =
       Goal;

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

(* eqn (5) *)
meyer_feta := y -> 1/2*(1 + (1 - y^2)*log((1 + y)/abs(1 - y))/(2*y)):

(* eqn (7) *)
meyer_lambda := y -> (1 - meyer_feta(y)) / (3 * y^2 * meyer_feta(y)):

(* enhancement factor from eqn (1) *)
meyer_f := x -> 1 + meyer_lambda(X2S*x/6)*x^2/(8*K_FACTOR_C):

f := (rs, zeta, xt, xs0, xs1) -> gga_kinetic(meyer_f, rs, zeta, xs0, xs1):


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_gga_x *)

c1 := 0.003:

f := x -> 1 + c1/X_FACTOR_C * x^2:

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
make_platonic := proc()
 global platonic_complexes,tetrahedron_complex,cube_complex,
   octahedron_complex,dodecahedron_complex,icosahedron_complex;
 local AC,TC,CC,OC,DC,IC,t0,t1,t2,s0,s1,AV,p,q,u,v,ui,s,t,n,m,d,M,A5,dp,nm,dd,r,k,i,E,g,
       E0,E1,E2,F1,F2,F3,C0,P0,PP0,PM,pole_colours,c;

 AC := [seq(table(),i=1..5)];
 IC := eval(AC[1]); OC := eval(AC[2]); TC := eval(AC[3]); CC := eval(AC[4]); DC := eval(AC[5]);

 t0 := (sqrt(5) + 1)/2;
 t1 := (sqrt(5) - 1)/2;
 t2 := expand(t1^2);
 s0 := sqrt(2/(5 + sqrt(5)));
 s1 := sqrt(2/(5 - sqrt(5)));

 dp := (u,v) -> simplify(expand(add(u[i] * v[i],i=1..3)));
 nm := (u) -> simplify(sqrt(dp(u,u)));
 dd := (u,v) -> nm(u -~ v);

 AV := [
 [
  [  0, t2, t1],[  0, t2,-t1],[  0,-t2, t1],[  0,-t2,-t1],
  [ t1,  0, t2],[-t1,  0, t2],[ t1,  0,-t2],[-t1,  0,-t2],
  [ t2, t1,  0],[ t2,-t1,  0],[-t2, t1,  0],[-t2,-t1,  0]
 ],
 [
  [ 1, 0, 0],[ 0, 1, 0],[ 0, 0, 1],
  [ 0, 0,-1],[ 0,-1, 0],[-1, 0, 0]
 ],
 [
  [ 1, 1, 1],[ 1,-1,-1],[-1, 1,-1],[-1,-1, 1]
 ],
 [
  [ 1, 1, 1],[ 1,-1,-1],[-1, 1,-1],[-1,-1, 1],
  [-1, 1, 1],[ 1,-1, 1],[ 1, 1,-1],[-1,-1,-1]
 ],
 [
  [  1,  1,  1],[  1, -1, -1],[ -1,  1, -1],[ -1, -1,  1],
  [  0, t0, t1],[  0, t0,-t1],[  0,-t0, t1],[  0,-t0,-t1],
  [ t1,  0, t0],[-t1,  0, t0],[ t1,  0,-t0],[-t1,  0,-t0],
  [ t0, t1,  0],[ t0,-t1,  0],[-t0, t1,  0],[-t0,-t1,  0],
  [ -1,  1,  1],[  1, -1,  1],[  1,  1, -1],[ -1, -1, -1]
 ]];

 for p from 1 to 5 do
  n := nops(AV[p]);
  AC[p]["num_vertices"] := n;
  AC[p]["vertices"] := [seq(i,i=1..AC[p]["num_vertices"])];
  r := simplify(sqrt(dp(AV[p][1],AV[p][1])));
  AC[p]["vertex_radius"] := r;
  AC[p]["embedding"] := table([seq(i = AV[p][i],i=1..n)]);
  AC[p]["normalised_embedding"] := table([seq(i = evalf(AV[p][i] /~ r),i=1..n)]);
  AC[p]["embedding_index"] := table():
  for i from 1 to n do
   AC[p]["embedding_index"][AV[p][i]] := i;
  od:
 od:

 for p from 1 to 5 do
  q := 6 - p;
  u := AV[p]; v := AV[q]; n := nops(u); m := nops(v);
  k := `if`(p = 5,"pentagons",`if`(p = 4,"squares","faces"));
  AC[p][k] := [seq(select(i -> dp(u[i],v[j]) = -1,[seq(i,i=1..n)]),j = 1..m)];
  AC[p][cat("oriented_",k)] :=
   [seq(orient_face(f,AC[p]["embedding"]),f in AC[p][k])];
  AC[p]["facets"] := AC[p][k];
  AC[p]["oriented_facets"] := AC[p][cat("oriented_",k)];
  AC[p][cat("num_",k)] := m;
  AC[p]["num_facets"] := m;
  AC[p]["plot"] := display(
   seq(polygon(map(j -> u[j],AC[p]["oriented_facets"][i]),colour=standard_colour(i)),i=1..m),
   scaling=constrained,axes=none);
  AC[p]["grey_plot"] := display(
   seq(polygon(map(j -> u[j],AC[p]["oriented_facets"][i]),colour=grey),i=1..m),
   scaling=constrained,axes=none);
  AC[p]["wireframe_plot"] := display(
   seq(polygon(map(j -> u[j],AC[p]["oriented_facets"][i]),colour=grey,style=wireframe),i=1..m),
   scaling=constrained,axes=none);
 od;

 for p from 1 to 3 do
  u := AV[p];
  E := {op(map(op,map(f -> [[f[1],f[2]],[f[1],f[3]],[f[2],f[3]]],AC[p]["faces"])))};
  AC[p]["edges"] := sort([op(E)]);
  AC[p]["short_edges"] := AC[p]["edges"];
  AC[p]["long_edges"] := [];
  r := dd(u[E[1][1]],u[E[1][2]]);
  AC[p]["short_edge_length"] := r;
 od:

 u := AV[4];
 r := 2;
 CC["short_edge_length"] := r;
 CC["short_edges"] :=
  select(e -> dd(u[e[1]],u[e[2]]) = r,[seq(seq([i,j],j=i+1..8),i=1..8)]);
 CC["long_edges"] := [[1,2],[1,3],[1,4],[5,8],[6,8],[7,8]];
 CC["edges"] := sort([op(CC["short_edges"]),op(CC["long_edges"])]);
 CC["faces"] := [];

 u := AV[5];
 r := sqrt(5) - 1;
 DC["short_edge_length"] := r;
 DC["short_edges"] :=
  select(e -> dd(u[e[1]],u[e[2]]) = 2,[seq(seq([i,j],j=i+1..20),i=1..20)]);

 A5 := eval(alternating_five);

 for p from 1 to 5 do
  u := AV[p];
  ui := AC[p]["embedding_index"];
  AC[p]["rotation_group"] := NULL;
  AC[p]["vertex_permutation"] := table();
  for g in A5["elements"] do
   M := Matrix(A5["matrix"][g]);
   s := [seq(ui[expand(convert(M . Vector(u[i]),list))],i=1..nops(u))];
   if type(s,list(posint)) then
    AC[p]["rotation_group"] := AC[p]["rotation_group"],g;
    AC[p]["vertex_permutation"][g] := s;
   fi;
  od:
  AC[p]["rotation_group"] := [AC[p]["rotation_group"]];
 od:

 for p from 2 to 4 do
  AC[p]["poles"] := table([2 = {}, 3 = {}, 5 = {}]);

  for g in AC[p]["rotation_group"] do
   if g = [1,2,3,4,5] then continue; fi;
   d := A5["element_order"][g];
   AC[p]["poles"][d] := {op(AC[p]["poles"][d]),A5["axis"][g],-~ A5["axis"][g]};
  od;

  for d in [2,3,5] do
   AC[p]["poles"][d] := [op(AC[p]["poles"][d])];
  od:
 od:

 PM := [1,-1];
 IC["poles"] := table():
 IC["poles"][2] := [
  [1,0,0],[0,1,0],[0,0,1],[-1,0,0],[0,-1,0],[0,0,-1],
  seq(seq(seq([i,j * t1,k * t0]/~2,k in PM),j in PM),i in PM),
  seq(seq(seq([k * t0,i,j * t1]/~2,k in PM),j in PM),i in PM),
  seq(seq(seq([j * t1,k * t0,i]/~2,k in PM),j in PM),i in PM)
 ];

 IC["poles"][3] := expand([seq(DC["embedding"][i] *~ sqrt(3)/3,i=1..20)]);

 IC["poles"][5] := [
  [  0, s0, s1],[  0, s0,-s1],[  0,-s0, s1],[  0,-s0,-s1],
  [ s1,  0, s0],[-s1,  0, s0],[ s1,  0,-s0],[-s1,  0,-s0],
  [ s0, s1,  0],[ s0,-s1,  0],[-s0, s1,  0],[-s0,-s1,  0]
 ];

 DC["poles"] := copy(IC["poles"]);

 E0 := NULL:
 E1 := NULL:
 E2 := NULL:
 F1 := NULL:
 F2 := NULL:
 F3 := NULL:
 for s in DC["rotation_group"] do
  t := DC["vertex_permutation"][s];
  E0 := E0,sort([t[1],t[ 5]]);
 od:
 for s in CC["rotation_group"] do
  t := DC["vertex_permutation"][s];
  E1 := E1,sort([t[1],t[10]]);
  E2 := E2,sort([t[1],t[17]]);
  F1 := F1,sort([t[1],t[ 9],t[10]]);
  F2 := F2,sort([t[1],t[ 5],t[17]]);
  F3 := F3,sort([t[1],t[10],t[17]]);
 od:
 E0 := sort([op({E0})]);
 E1 := sort([op({E1})]);
 E2 := sort([op({E2})]);
 F1 := sort([op({F1})]);
 F2 := sort([op({F2})]);
 F3 := sort([op({F3})]);

 DC["short_edges"] := E0;
 DC["long_edges_a"] := E1;
 DC["long_edges_b"] := E2;
 DC["edges"] := sort([op(E0),op(E1),op(E2)]);

 DC["faces_a"] := F1;
 DC["faces_b"] := F2;
 DC["faces_c"] := F3;
 DC["faces"] := sort([op(F1),op(F2),op(F3)]);

 DC["inscribed_cubes"] := [
  [3, 5, 8, 10, 11, 13, 16, 18],
  [4, 5, 8,  9, 12, 14, 15, 19],
  [2, 6, 7,  9, 12, 13, 16, 17],
  [1, 6, 7, 10, 11, 14, 15, 20],
  [1, 2, 3,  4, 17, 18, 19, 20]
 ];

 v := DC["embedding"];

 DC["all_edges_plot"] :=
  display(
   seq(line(v[e[1]],v[e[2]],colour=blue) ,e in DC["short_edges"]),
   seq(line(v[e[1]],v[e[2]],colour=red)  ,e in DC["long_edges_a"]),
   seq(line(v[e[1]],v[e[2]],colour=green),e in DC["long_edges_b"]),
   scaling=constrained,axes=none
  );

 DC["triangle_plot"] :=
  display(
   seq(polygon([v[e[1]],v[e[2]],v[e[3]]],colour=blue ),e in DC["faces_a"]),
   seq(polygon([v[e[1]],v[e[2]],v[e[3]]],colour=red  ),e in DC["faces_b"]),
   seq(polygon([v[e[1]],v[e[2]],v[e[3]]],colour=green),e in DC["faces_c"]),
   scaling=constrained,axes=none
  );

 PP0 := NULL:
 for k from 1 to 5 do
  C0 := DC["inscribed_cubes"][k];
  E0 := [seq(seq([C0[i],C0[j]],j=i+1..8),i=1..7)]:
  E0 := select(e -> simplify(expand(dd(v[e[1]],v[e[2]]))) = 2,E0);
  P0 := display(seq(line(v[e[1]],v[e[2]],colour=standard_colour(k)),e in E0),
                scaling=constrained,axes=none);
  PP0 := PP0,P0;
 od:
 DC["cubes_plot"] := display(DC["grey_plot"],PP0);

 pole_colours := table([2 = green, 3 = red, 5 = blue]):

 IC["pole_length"] := table([2 = 0.8, 3 = 0.8, 5 = 0.9]):
 OC["pole_length"] := table([2 = 1.2, 3 = 0.8, 5 = 0.9]):
 TC["pole_length"] := table([2 = 1.3, 3 = 1.2, 5 = 0.9]):
 CC["pole_length"] := table([2 = 1.3, 3 = 2.0, 5 = 0.9]):
 DC["pole_length"] := table([2 = 2.0, 3 = 2.0, 5 = 1.7]):

 for p from 1 to 5 do
  AC[p]["unit_pole_plot"] := table():
  AC[p]["pole_plot"] := table():
  for d in [2,3,5] do
   r := AC[p]["pole_length"][d];
   c := pole_colours[d];
   AC[p]["pole_plot"][d] :=
     map(u -> line(-r *~ u,r *~u,colour = c,thickness=4),AC[p]["poles"][d]);
   AC[p]["unit_pole_plot"][d] :=
     map(u -> line(-1.1 *~ u,1.1 *~u,colour = c,thickness=4),AC[p]["poles"][d]);
  od:
  AC[p]["all_poles_plot"] :=
   display(seq(op(AC[p]["pole_plot"][d]),d in [2,3,5]),
	   scaling=constrained,axes=none);
  AC[p]["all_unit_poles_plot"] :=
   display(seq(op(AC[p]["unit_pole_plot"][d]),d in [2,3,5]),
	   scaling=constrained,axes=none);
  if p >= 2 and p <= 4 then
   AC[p]["sample_poles_plot"] :=
    display(AC[p]["pole_plot"][2][1],
	    AC[p]["pole_plot"][3][1],
	    scaling=constrained,axes=none);
  else
   AC[p]["sample_poles_plot"] :=
    display(AC[p]["pole_plot"][2][1],
	    AC[p]["pole_plot"][3][1],
	    AC[p]["pole_plot"][5][1],
	    scaling=constrained,axes=none);
  fi:
 od:

 IC["dual_factor"] := 3;
 OC["dual_factor"] := 3;
 TC["dual_factor"] := 3;
 CC["dual_factor"] := 1;
 DC["dual_factor"] := 5 - 2 * sqrt(5);

 for p from 1 to 5 do
  q := 6 - p;
  v := AC[q]["embedding"];
  m := nops(AC[q]["oriented_facets"]);
  r := AC[p]["dual_factor"];

  AC[p]["dual_plot"] := display(
   AC[p]["wireframe_plot"],
   seq(polygon(map(j -> evalf(v[j] /~ (-r)),AC[q]["oriented_facets"][i]),colour=standard_colour(i)),i=1..m),
   scaling=constrained,axes=none
  );
 od:

 platonic_complexes   := eval(AC);
 icosahedron_complex  := eval(IC);
 octahedron_complex   := eval(OC);
 tetrahedron_complex  := eval(TC);
 cube_complex         := eval(CC);
 dodecahedron_complex := eval(DC);
 NULL;
end:


text:
program ifst;
var ch : char;
begin
    readln(ch);
    if ch = 'a' then
        writeln('It is ''a'' ')
    else
        writeln('It is not ''a'' ')
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


text:
//impulse: /give @p spawn_egg 1 0 {display:{Name:"Cleansing Gel Spawner"},EntityTag:{id:"minecraft:chicken",CustomName:"MMH_CleansingGelSpawner",Silent:1,NoGravity:1}}
//MMH_loadChunks()
///clone ~1 ~1 ~1 ~3 ~5 ~3 24 1 24

#MMH
repeat process MMH_cleansingGelSpawner {
if: /testfor @e[type=Chicken,name=MMH_CleansingGelSpawner]
then {
  MMH_loadChunks()
  /execute @e[type=Chicken,name=MMH_CleansingGelSpawner] ~ ~ ~ clone 24 1 24 26 5 26 ~-1 ~ ~-1 masked
  /execute @e[type=Chicken,name=MMH_CleansingGelSpawner] ~ ~ ~ fill ~ ~ ~ ~ ~ ~ air 0 replace pumpkin
  /tp @e[type=Chicken,name=MMH_CleansingGelSpawner] ~ -100 ~
}
}

///summon item ~-1 ~-1 ~-1 {CustomName:"ACV_CleansingGel",PickupDelay:32767,Item:{Count:1,Damage:1,id:"minecraft:wool"}}


text:
`delta/F3` := (N::posint) -> (x) -> [
 `delta/F`(N)({1,2,3})(1,2)(x),
 `delta/F`(N)({1,2,3})(2,3)(x),
 `delta/F`(N)({1,2,3})(3,1)(x)
];

`Delta/F3` := (N::posint) -> (u) -> [u,u];

`f_plus/F3` := (N::posint) -> proc(vw)
 local v,w;
 v,w := op(vw);
 return [-~ v,v,v -~ w];
end:

`f_minus/F3` := (N::posint) -> proc(vw)
 local v,w;
 v,w := op(vw);
 return [-~ v,v,w -~ v];
end:

`g_plus/F3`  := (N::posint) -> proc(vw)
 local v,w;
 v,w := op(vw);
 return [v,-~ w, -~ v];
end:

`g_minus/F3`  := (N::posint) -> proc(vw)
 local v,w;
 v,w := op(vw);
 return [v,-~ v, -~ w];
end:

######################################################################

`is_element/X/F3` := (N::posint) -> proc(xyz)
 global reason;
 local tag,x,y,z;

 tag := "is_element/X/F3";

 if not(type(xyz,list) and nops(xyz) = 3) then
  reason := [tag,"xyz is not a list of length 3",xyz];
  return false;
 fi;

 x,y,z := op(xyz);
 if not(`is_element/R`(N)(x) and
        `is_element/R`(N)(y) and
        `is_element/R`(N)(z)) then
  reason := [tag,"x, y and z are not vectors in R^N",x,y,z];
  return false;
 fi;

 if (x -~ y = [0$N] or y -~ z = [0$N] or z -~ x = [0$N]) then
  reason := [tag,"x, y and z are not distinct",x,y,z];
  return false;
 fi;

 if x +~ y <> [0$N] then
  reason := [tag,"x+y <> 0",x,y];
  return false;
 fi;

 if `norm_2/R`(N)(x) <> 1 or `norm_2/R`(N)(y) <> 1 then
  reason := [tag,"x and y are not normalised",x,y];
  return false;
 fi;

 return true;
end:

`is_equal/X/F3` := (N::posint) -> (xyz0,xyz1) ->
  evalb(simplify(xyz0 -~ xyz1) = [[0$N]$3]);

`is_leq/X/F3` := NULL;
`list_elements/X/F3` := NULL;
`count_elements/X/F3` := NULL;

`random_element/X/F3` := (N::posint) -> proc()
 local x,y,z;
 x := `random_element/sphere`(N-1)();
 y := -~ x;
 z := `random_element/R`(N)();
 while z = x or z = y do
  z := `random_element/R`(N)();
 od;
 return [x,y,z];
end:

######################################################################

`is_element/Y/F3` := (N::posint) -> proc(xyz)
 local tag,x,y,z;
 global reason;

 tag := "is_element/Y/F3";

 if not(`is_element/X/F3`(N)(xyz)) then
  reason := [tag,"xyz in not in X",xyz,reason];
  return false;
 fi;

 x,y,z := op(xyz);
 if `d_2/R`(N)(z,x) <> 1 and `d_2/R`(N)(z,y) <> 1 then
  reason := [tag,"z is not at distance 1 from x or y",x,y,z];
  return false;
 fi;

 return true;
end:

`is_equal/Y/F3` := (N::posint) -> (xyz0,xyz1) ->
  evalb(simplify(xyz0 -~ xyz1) = [[0$N]$3]);

`is_leq/Y/F3` := NULL;
`list_elements/Y/F3` := NULL;
`count_elements/Y/F3` := NULL;

`random_element/Y/F3` := (N::posint) -> proc()
 local x,y,w,z;
 x := `random_element/sphere`(N-1)();
 y := -~ x;
 w := `random_element/sphere`(N-1)();
 if rand(0..1)() = 0 then
  z := x +~ w;
 else
  z := y +~ w;
 fi;
 return [x,y,z];
end:

######################################################################

`is_element/D_plus/F3` := (N::posint) -> proc(xyz)
 local tag,x,y,z;
 global reason;
 tag := "is_element/D_plus/F3";

 if not(`is_element/X/F3`(N)(xyz)) then
  reason := [tag,"xyz in not in X",xyz,reason];
  return false;
 fi;

 x,y,z := op(xyz);
 if `d_2/R`(N)(z,y) > 1 then
  reason := [tag,"z is not at distance <= 1 from y",x,y,z];
  return false;
 fi;

 return true;
end:

`is_equal/D_plus/F3` := (N::posint) -> (xyz0,xyz1) ->
  evalb(simplify(xyz0 -~ xyz1) = [[0$N]$3]);

`is_leq/D_plus/F3` := NULL;
`list_elements/D_plus/F3` := NULL;
`count_elements/D_plus/F3` := NULL;

`random_element/D_plus/F3` := (N::posint) -> proc()
 local x,y,w,z;
 x := `random_element/sphere`(N-1)();
 y := -~ x;

 if rand(0..5)() = 0 then
  z := [0$N];
 else
  w := `random_element/sphere`(N-1)() *~ rand(1..5)()/5;
  z := y +~ w;
 fi;
 return [x,y,z];
end:

######################################################################

`is_element/D_minus/F3` := (N::posint) -> proc(xyz)
 local tag,x,y,z;
 global reason;

 tag := "is_element/D_minus/F3";

 if not(`is_element/X/F3`(N)(xyz)) then
  reason := [tag,"xyz in not in X",xyz,reason];
  return false;
 fi;

 x,y,z := op(xyz);
 if `d_2/R`(N)(z,x) > 1 then
  reason := [tag,"z is not at distance <= 1 from x",x,y,z];
  return false;
 fi;

 return true;
end:

`is_equal/D_minus/F3` := (N::posint) -> (xyz0,xyz1) ->
  evalb(simplify(xyz0 -~ xyz1) = [[0$N]$3]);

`is_leq/D_minus/F3` := NULL;
`list_elements/D_minus/F3` := NULL;
`count_elements/D_minus/F3` := NULL;

`random_element/D_minus/F3` := (N::posint) -> proc()
 local x,y,w,z;
 x := `random_element/sphere`(N-1)();
 y := -~ x;

 if rand(0..5)() = 0 then
  z := [0$N];
 else
  w := `random_element/sphere`(N-1)() *~ rand(1..5)()/5;
  z := x +~ w;
 fi;
 return [x,y,z];
end:

######################################################################

`is_element/D_zero/F3` := (N::posint) -> proc(xyz)
 local tag,x,y,z;
 global reason;

 tag := "is_element/D_zero/F3";

 if not(`is_element/X/F3`(N)(xyz)) then
  reason := [tag,"xyz in not in X",xyz,reason];
  return false;
 fi;

 x,y,z := op(xyz);
 if `d_2/R`(N)(z,x)^2 < 1 or `d_2/R`(N)(z,y)^2 < 1 then
  reason := [tag,"z is not at distance >= 1 from x and y",x,y,z];
  return false;
 fi;

 return true;
end:

`is_equal/D_zero/F3` := (N::posint) -> (xyz0,xyz1) ->
  evalb(simplify(xyz0 -~ xyz1) = [[0$N]$3]);

`is_leq/D_zero/F3` := NULL;
`list_elements/D_zero/F3` := NULL;
`count_elements/D_zero/F3` := NULL;

`random_element/D_zero/F3` := (N::posint) -> proc()
 local x,y,w,z,r,i;
 x := `random_element/sphere`(N-1)();
 y := -~ x;

 z := x;
 r := 0;
 while r < 1 do
  if rand(0..5)() = 0 then
   w := `random_element/sphere`(N-1)();
   if rand(0..1)() = 0 then
    z := x +~ w;
   else
    z := y +~ w;
   fi;
  else
   z := `random_element/sphere`(N-1)(100) *~ (rand(1..300)()/100);
  fi;
  r := min(add((x[i]-z[i])^2,i=1..N),add((y[i]-z[i])^2,i=1..N));
 od:

 return [x,y,z];
end:

######################################################################

`q_plus/F3` := (N::posint) -> proc(xyz)
 local x,y,z;
 x,y,z := op(xyz);
 return y +~ (z -~ y) /~ `norm_2/R`(N)(z -~ y);
end:

######################################################################

`q_minus/F3` := (N::posint) -> proc(xyz)
 local x,y,z;
 x,y,z := op(xyz);
 return y +~ (z -~ x) /~ `norm_2/R`(N)(z -~ x);
end:

######################################################################

`q_zero/F3` := (N::posint) -> proc(xyz)
 local x,y,z;
 x,y,z := op(xyz);
 if z = [0$N] then
  return [0$N];
 else
  return (2*`dot/R`(N)(y,z)/`norm_2/R`(N)(z)^2) *~ z;
 fi;
end:

######################################################################

`q/F3` := (N::posint) -> proc(xyz)
 local x,y,z;
 x,y,z := op(xyz);
 if   `is_element/D_plus/F3`(N)(xyz) then
  return `q_plus/F3`(N)(xyz);
 elif `is_element/D_minus/F3`(N)(xyz) then
  return `q_minus/F3`(N)(xyz);
 else
  return `q_zero/F3`(N)(xyz);
 fi;
end:

`r/F3` := (N::posint) -> proc(xyz)
 local x,y,z,x0,y0,z0,n0;
 x,y,z := op(xyz);

 x0 := (x -~ y)/~2;
 n0 := `norm_2/R`(N)(x0);
 x0 := x0 /~ n0;
 y0 := -~ x0;
 z0 := (z -~ ((x +~ y)/~2))/~ n0;
 return [x0,y0,`q/F3`(N)([x0,y0,z0])];
end:


text:
######################################################################

`is_element/ord_simplex_interior` := (A::set) -> proc(Rx)
 local R,x;
 global reason;

 if not(type(Rx,list) and nops(Rx) = 2) then
  reason := [convert(procname,string),"Rx cannot be split as [R,x]",Rx];
  return false;
 fi;

 R,x := op(Rx);

 if not(`is_element/ord`(A)(R)) then
  reason := [convert(procname,string),"R is not an order on A",R,A,reason];
  return false;
 fi;

 if not(`is_element/simplex_interior`(A)(x)) then
  reason := [convert(procname,string),"x is not in the interior of the A-simplex",x,A,reason];
  return false;
 fi;

 return true;
end;

`is_equal/ord_simplex_interior` := (A::set) -> proc(Rx,Sy)
 global reason;

 if not `is_equal/ord`(A)(Rx[1],Sy[1]) then
  reason := [convert(procname,string),"R <> S",Rx[1],Sy[1],reason];
  return false;
 fi;

 if not `is_equal/simplex`(A)(Rx[2],Sy[2]) then
  reason := [convert(procname,string),"x <> y",Rx[2],Sy[2],reason];
  return false;
 fi;

 return true;
end;

`is_leq/ord_simplex_interior` := NULL;

`random_element/ord_simplex_interior` := (A::set) -> proc(d::posint := 5)
 [`random_element/ord`(A)(),
  `random_element/simplex_interior`(A)(d)];
end;

`list_elements/ord_simplex_interior` := NULL;
`count_elements/ord_simplex_interior` := NULL;

`phi/ord_simplex_interior/one_cubes_prime` := (A::set) -> proc(Rx)
 local R,x,y,n,i,j,p,q;

 R,x := op(Rx);
 y := table();
 n := nops(R);

 for i from 1 to n do
  p := add(x[R[j]],j=1..i-1);
  q := p + x[R[i]];
  y[R[i]] := [[p],[q]];
 od;

 return [R,eval(y)];
end;


text:
(* type: lda_exc *)

a :=  0.0311:
b := -0.048:
c :=  0.009:
d := -0.017:

f_ls := (rs, zeta) -> a*log(rs) + b + c*rs*log(rs) + d*rs:
f_ms := (rs, zeta) -> -0.06156 + 0.01898*log(rs):
f_hs := (rs, zeta) -> 0.438/rs + 1.325/rs^(3/2) - 1.47/rs^2 - 0.4/rs^(5/2):

f := (rs, zeta) -> my_piecewise3(
  rs < 0.7, f_ls(rs, zeta),  my_piecewise3(rs < 10,  f_ms(rs, zeta), f_hs(rs, zeta))):


text:
program rcall;
procedure f; begin call f end;
begin call f end.

max_issues_repo_name: shouth/LanguageProcessing


text:
`is_element/permutations` := (n::nonnegint) -> proc(s)
 type(s,list(posint)) and
 nops(s) = n and
 min(op(s)) = 1 and
 max(op(s)) = n and
 nops({op(s)}) = n;
end:

`is_equal/permutations` := (n) -> (s,t) -> evalb(s = t);

`is_leq/permutations` := NULL;

`random_element/permutations` := (n::nonnegint) ->
 combinat[randperm](n);

`list_elements/permutations` := (n::nonnegint) ->
 combinat[permute](n);

`count_elements/permutations` := (n::nonnegint) -> n!;

`id/permutations` := proc(n::nonnegint) local i; [seq(i,i=1..n)]; end:

`o/permutations` := (n::nonnegint) -> proc()
 local i;
 apply_assoc((s,t) -> [seq(s[t[i]],i=1..n)],
             [seq(i,i=1..n)])(args);
end:

`inv/permutations` := (n::nonnegint) -> proc(s)
 local i,t;
 t := table();
 for i from 1 to n do t[s[i]] := i; od;
 return [seq(t[i],i=1..n)];
end:

`length_set/permutations` := (n::nonnegint) -> proc(s)
 local L,i,j;

 L := [seq(seq([i,j],j=i+1..n),i=1..n-1)];
 L := select(ij -> s[ij[1]] > s[ij[2]],L);
 return L;
end:

`unordered_length_set/permutations` := (n::nonnegint) -> (s) ->
  map(u -> {op(u)},{op(`length_set/permutations`(n)(s))});

`length/permutations` := (n::nonnegint) -> proc(s)
 nops(`length_set/permutations`(n)(s));
end:

`sgn/permutations` := (n::nonnegint) -> proc(s)
 local i,j;
 return signum(mul(mul(s[j]-s[i],j=i+1..n),i=1..n-1));
end:

`is_even/permutations` := (n::nonnegint) -> (s) ->
 evalb(`sgn/permutations`(n)(s) =  1);

`is_odd/permutations` := (n::nonnegint) -> (s) ->
 evalb(`sgn/permutations`(n)(s) = -1);

`switch/permutations` := (n::nonnegint) -> proc(i::posint)
  local j;
  if i >= n then error("Should have 1 <= i < n"); fi;

  return [seq(j,j=1..i-1),i+1,i,seq(j,j=i+2..n)];
end:

`t/permutations` := (n::nonnegint) -> proc(k::posint)
 local i;
 [seq(i,i=1..k-1),seq(i+1,i=k..n-1),k];
end:

`t_inv/permutations` := (n::nonnegint) -> proc(k::posint)
 local i;
 [seq(i,i=1..k-1),n,seq(i-1,i=k+1..n)];
end:

`t_word/permutations` := (n::nonnegint) -> proc(k::posint)
 local i;
 [seq(i,i=k..n-1)];
end:

`to_coxeter_word/permutations` := (n::nonnegint) -> proc(w)
 local m,v;

 if n <= 1 then return []; fi;

 m := w[n];
 v := `o/permutations`(n)(`t_inv/permutations`(n)(m),w);
 v := [op(1..(n-1),v)];
 return [op(`t_word/permutations`(n)(m)),
         op(`to_coxeter_word/permutations`(n-1)(v))];
end:

`from_coxeter_word/permutations` := (n::nonnegint) -> proc(x)
 `o/permutations`(n)(op(map(`switch/permutations`(n),x)));
end:

# Here are two kinds of Coxeter words that reduce to the identity.
# They are used in the proof of the completeness of the Coxeter relations.

`long_coxeter_rel0/permutations` := (n) -> proc(l)
 local i;
 [seq(i,i=l..n-1),seq(i,i=l..n-1),seq(i,i=n-2..l,-1),seq(i,i=n-1..l+1,-1)];
end:

`long_coxeter_rel1/permutations` := (n) -> proc(k,l)
 local i;
 `if`(k <= l,
       [seq(i,i=k..n-1),seq(i,i=l..n-1),seq(i,i=n-2..k,-1),seq(i,i=n-1..l+1,-1)],
       [seq(i,i=k..n-1),seq(i,i=l..n-1),seq(i,i=n-2..k-1,-1),seq(i,i=n-1..l,-1)]);
end:

`coxeter_reduce_once/permutations` := (n) -> proc(x)
 local i,j,k,l,l0,l1,m,p,d,y,z,i0,p0,d0,z0,z1,z2;

 # Move s[p] left of s[q] if p < q - 1
 y := x;
 m := nops(y);
 z := [];
 i := 1;
 while i <= m do
  j := i + 1;
  while j <= m and y[j] < y[i] - 1 do
   z := [op(z),y[j]];
   j := j + 1;
  od;
  z := [op(z),y[i]];
  i := j;
 od;
 y := z;
 m := nops(y);

 z := [];
 i := 1;
 j := 1;
 while i <= m do
  while (j < m and y[j+1] = y[i]) do j := j + 1; od;
  if modp(j - i + 1,2) = 1 then z := [op(z),y[i]]; fi;
  i := j + 1;
  j := i;
 od;
 y := z;
 m := nops(y);

 d := table():
 l := table():
 i0 := 0;
 l0 := 0;

 for i from 1 to m do
  if i < m and abs(y[i+1] - y[i]) = 1 then
   d[i] := y[i+1] - y[i];
   j := i + 1;
   while j <= m and y[j] = y[i] + modp(j - i,2) * d[i] do j := j + 1; od;
   l[i] := j - i;
  else
   d[i] := 0;
   l[i] := 1;
  fi;
  if l[i] > l0 then
   i0 := i;
   l0 := l[i];
  fi;
 od:

 if l0 < 3 then return y; fi;

 z0 := [seq(y[i],i=1..i0-1)];
 z2 := [seq(y[i],i=i0+l0..m)];
 l1 := mods(l0,6);
 p0 := y[i0];
 d0 := d[i0];

 if   l1 = 0             then z1 := [];
 elif l1 = 1             then z1 := [p0];
 elif l1 = 2             then z1 := [p0,p0+d0];
 elif l1 = 3 and d0 =  1 then z1 := [p0,p0+1,p0];
 elif l1 = 3 and d0 = -1 then z1 := [p0-1,p0,p0-1];
 elif l1 = -2            then z1 := [p0+d0,p0];
 elif l1 = -1            then z1 := [p0+d0];
 fi;

 return [op(z0),op(z1),op(z2)];
end:

`coxeter_reduce/permutations` := (n) -> proc(w)
 local u,v;
 u := NULL;
 v := w;
 while u <> v do
  u := v;
  v := `coxeter_reduce_once/permutations`(n)(u);
 od:
 return u;
end:


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_mgga_x *)

$include "mgga_x_gx.mpl"

mmu := 0.001015549:
pbe_gx := x -> 1/(1 + mmu*x^2):

f := (rs, x, t, u) ->
  f_gx_a(malpha(x, t)) * pbe_gx(x):


text:
flavor 1
srclang 3
id 65535
numfuncs 4
import "Arith.mplt"
import "/home/bravewtz/Desktop/openarkcompiler/libjava-core/libjava-core.mplt"
entryfunc &LArith_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V
fileinfo {
  @INFO_filename "Arith.jar"}
srcfileinfo {
  1 "Arith.java"}
type $__class_meta__ <struct {
  @shadow <* void> final,
  @monitor i32 final,
  @classloader u16 final,
  @objsize u16 final,
  @itab <* void> final,
  @vtab <* void> final,
  @gctib <* void> final,
  @classinforo <* void> final,
  @clinitbridge <* void> final}>
type $__class_meta_ro__ <struct {
  @classname <* void> final,
  @ifields <* void> final,
  @methods <* void> final,
  @superclass_or_componentclass <* void> final,
  @numoffields u16 final,
  @numofmethods u16 final,
  @flag u16 final,
  @numofsuperclasses u16 final,
  @padding u32 final,
  @mod i32 final,
  @annotation i32 final,
  @clinitAddr i32 final}>
type $__method_info__ <struct {
  @method_in_vtab_index i64 final,
  @declaringclass i64 final,
  @addr i64 final,
  @mod i32 final,
  @methodname i32 final,
  @signaturename i32 final,
  @annotationvalue i32 final,
  @flag u16 final,
  @argsize u16 final,
  @padding u32 final}>
type $__method_info_compact__ <struct {
  @method_in_vtab_index i32 final,
  @addr i32 final,
  @lebPadding0 u8 final}>
type $__field_info__ <struct {
  @offset u64 final,
  @mod u32 final,
  @flag u16 final,
  @index u16 final,
  @typeName i64 final,
  @fieldname u32 final,
  @annotation u32 final,
  @declaringclass <* void> final}>
type $__field_info_compact__ <struct {
  @offset u32 final,
  @lebPadding0 u8 final}>
type $__superclass_meta__ <struct {
  @superclassinfo <* void> final}>
type $MUIDFuncDefTabEntry <struct {
  @funcUnifiedAddr u64}>
type $MUIDFuncInfTabEntry <struct {
  @funcSize u32,
  @funcName u32}>
type $MUIDFuncDefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDDataDefTabEntry <struct {
  @dataUnifiedAddr u64}>
type $MUIDDataDefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDUnifiedUndefTabEntry <struct {
  @globalAddress u64}>
type $MUIDUnifiedUndefMuidTabEntry <struct {
  @muidLow u64,
  @muidHigh u64}>
type $MUIDRangeTabEntry <struct {
  @tabBegin <* void>,
  @tabEnd <* void>}>
javaclass $LArith_3B <$LArith_3B> public
var $Ljava_2Flang_2FSystem_3B_7Cout extern <* <$Ljava_2Fio_2FPrintStream_3B>> final public static
var $__cinf_Ljava_2Flang_2FString_3B <$__class_meta__>
func &MCC_GetOrInsertLiteral () <* <$Ljava_2Flang_2FString_3B>>
var $__vtb_LArith_3B fstatic <[11] <* void>> = [16, 24, 32, 36, 8, 4, 52, 12, 20, 48, 28]
var $__cinf_LArith_3B <$__class_meta__> public
var $__methods_info__LArith_3B fstatic <[4] <$__method_info__>> public = [[1= 0xfff6, 2= addrof ptr $__cinf_LArith_3B, 3= addroffunc ptr &LArith_3B_7Craise__sigfpe_7C_28_29V, 4= 266, 5= 4, 6= 56, 7= 72, 8= 0x2f41, 9= 0, 10= 0], [1= 0xfff6, 2= addrof ptr $__cinf_LArith_3B, 3= addroffunc ptr &LArith_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V, 4= 9, 5= 88, 6= 108, 7= 72, 8= 0xa201, 9= 1, 10= 0], [1= 0xfff6, 2= addrof ptr $__cinf_LArith_3B, 3= addroffunc ptr &LArith_3B_7C_3Cinit_3E_7C_28_29V, 4= 0x10001, 5= 200, 6= 56, 7= 72, 8= 0xad81, 9= 1, 10= 0], [1= 0xfff6, 2= addrof ptr $__cinf_LArith_3B, 3= addroffunc ptr &LArith_3B_7CTestMain_7C_28I_29I, 4= 9, 5= 228, 6= 264, 7= 72, 8= 0xc081, 9= 1, 10= 0]]
var $__cinf_Ljava_2Flang_2FObject_3B extern <$__class_meta__> public
var $__superclasses__LArith_3B fstatic <[1] <$__superclass_meta__>> public = [[1= 0x4000000000000000]]
var $__classinforo__LArith_3B fstatic <$__class_meta_ro__> public = [1= 7, 2= 0, 3= addrof ptr $__methods_info__LArith_3B, 4= addrof ptr $__superclasses__LArith_3B, 5= 0, 6= 4, 7= 0, 8= 1, 9= 0, 10= 33, 11= 72, 12= 0]
var $MCC_GCTIB__LArith_3B fstatic <* void> public
var $classStateInitialized u64
var $__cinf_LArith_3B <$__class_meta__> public = [1= 0x17eca191257f4b, 2= 0, 3= 0xffff, 4= 0, 5= 0, 6= addrof ptr $__vtb_LArith_3B, 7= addrof ptr $MCC_GCTIB__LArith_3B, 8= addrof ptr $__classinforo__LArith_3B, 9= addrof ptr $classStateInitialized]
var $__muid_classmetadata_bucket$$Arith_jar <[1] <* void>> public = [addrof ptr $__cinf_LArith_3B]
func &MCC_Reflect_ThrowCastException nosideeffect () void
func &MCC_Reflect_Check_Casting_NoArray nosideeffect () void
func &MCC_Reflect_Check_Casting_Array nosideeffect () void
func &MCC_CheckThrowPendingException nosideeffect () void
func &MCC_PreNativeCall (var %caller ref) <* void>
func &MCC_PostNativeCall (var %env <* void>) void
func &MCC_DecodeReference nosideeffect (var %obj ref) ref
func &MCC_CallFastNative (var %func <* void>) <* void>
func &MCC_CallSlowNative0 (var %func <* void>) <* void>
func &MCC_CallSlowNative1 (var %func <* void>) <* void>
func &MCC_CallSlowNative2 (var %func <* void>) <* void>
func &MCC_CallSlowNative3 (var %func <* void>) <* void>
func &MCC_CallSlowNative4 (var %func <* void>) <* void>
func &MCC_CallSlowNative5 (var %func <* void>) <* void>
func &MCC_CallSlowNative6 (var %func <* void>) <* void>
func &MCC_CallSlowNative7 (var %func <* void>) <* void>
func &MCC_CallSlowNative8 (var %func <* void>) <* void>
func &MCC_CallFastNativeExt (var %func <* void>) <* void>
func &MCC_CallSlowNativeExt (var %func <* void>) <* void>
func &MCC_SetReliableUnwindContext nosideeffect () void
var $__reg_jni_func_tab$$Arith_jar <[1] <* void>> = [-72057594037927936]
func &Java_Arith_raise_1sigfpe__ weak () void
func &MCC_CannotFindNativeMethod nosideeffect () <* void>
func &MCC_FindNativeMethodPtrWithoutException nosideeffect () <* void>
func &MCC_DummyNativeMethodPtr nosideeffect () <* void>
var $__reg_jni_tab$$Arith_jar <[1] <* void>> = [7, 39]
var $__cinf_Ljava_2Flang_2FSystem_3B extern <$__class_meta__>
func &MCC_getFuncPtrFromItabSecondHash64 nosideeffect () ptr
var $__cinf_Ljava_2Flang_2FArithmeticException_3B extern ptr
var $__muid_func_def_tab$$Arith_jar fstatic <[4] <$MUIDFuncDefTabEntry>> = [[1= addroffunc ptr &LArith_3B_7C_3Cinit_3E_7C_28_29V], [1= addroffunc ptr &LArith_3B_7Craise__sigfpe_7C_28_29V], [1= addroffunc ptr &LArith_3B_7CTestMain_7C_28I_29I], [1= addroffunc ptr &LArith_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V]]
var $__muid_func_inf_tab$$Arith_jar fstatic <[4] <$MUIDFuncInfTabEntry>> = [[1= addroffunc ptr &LArith_3B_7C_3Cinit_3E_7C_28_29V, 2= addroffunc ptr &LArith_3B_7C_3Cinit_3E_7C_28_29V], [1= addroffunc ptr &LArith_3B_7Craise__sigfpe_7C_28_29V, 2= addroffunc ptr &LArith_3B_7Craise__sigfpe_7C_28_29V], [1= addroffunc ptr &LArith_3B_7CTestMain_7C_28I_29I, 2= addroffunc ptr &LArith_3B_7CTestMain_7C_28I_29I], [1= addroffunc ptr &LArith_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V, 2= addroffunc ptr &LArith_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V]]
var $__muid_func_def_muid_tab$$Arith_jar fstatic <[4] <$MUIDFuncDefMuidTabEntry>> = [[1= -4001351761653078127, 2= -1310918331262075617], [1= -8478930883496526569, 2= -2136361627923192435], [1= 0x1ea56a56b0696218, 2= -2592354488807809467], [1= -4467004309115070999, 2= -2701549801737661972]]
var $__muid_func_muid_idx_tab$$Arith_jar fstatic <[4] u32> = [3, 2, 1, 0]
var $__muid_data_def_tab$$Arith_jar fstatic <[1] <$MUIDDataDefTabEntry>> = [[1= addrof ptr $__cinf_LArith_3B]]
var $__muid_data_def_muid_tab$$Arith_jar fstatic <[1] <$MUIDDataDefMuidTabEntry>> = [[1= -7613525522286655244, 2= -2842220972116417222]]
var $__muid_func_undef_tab$$Arith_jar fstatic <[13] <$MUIDUnifiedUndefTabEntry>> = [[1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0], [1= 0]]
var $__muid_func_undef_muid_tab$$Arith_jar fstatic <[13] <$MUIDUnifiedUndefMuidTabEntry>> = [[1= 0x3e32352aee789835, 2= -3887705395317205813], [1= -4187412136968710015, 2= -3800091941095621250], [1= 0x6742c234127e0a27, 2= -3762262047879347071], [1= 0x7230554331c55d92, 2= -3676689525926909155], [1= -715372855679083712, 2= -2647497990906227723], [1= -7464356948810446352, 2= -2259485500590180091], [1= 0x783627f2afd1cbde, 2= -2046851302095768916], [1= -2701934576591406938, 2= -1693831364093527548], [1= 0x126777a7fe39e1fb, 2= -1314856249532362766], [1= 0xbf40578f3343f7a, 2= -1198421541845410999], [1= 0xf83f176d111918a, 2= -976603133906039419], [1= 0x7ca2bdf69a6c7c94, 2= -801329978528900548], [1= 0x477aafa4d7dd102b, 2= -442561182569419835]]
var $__muid_data_undef_tab$$Arith_jar fstatic <[4] <$MUIDUnifiedUndefTabEntry>> = [[1= addrof ptr $__cinf_Ljava_2Flang_2FObject_3B], [1= addrof ptr $Ljava_2Flang_2FSystem_3B_7Cout], [1= addrof ptr $__cinf_Ljava_2Flang_2FArithmeticException_3B], [1= addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B]]
var $__muid_data_undef_muid_tab$$Arith_jar fstatic <[4] <$MUIDUnifiedUndefMuidTabEntry>> = [[1= -567417612161374449, 2= -3298852447504547670], [1= 0x191283ac418c4bb9, 2= -1676204161023949463], [1= -8918818566091393965, 2= -1601760931404013578], [1= -5921653145571052587, 2= -171150348656858163]]
var $__muid_range_tab$$Arith_jar fstatic <[29] <$MUIDRangeTabEntry>> = [[1= 0x33693ce2ebe87faa, 2= -3126084233764470422], [1= 0x46d1a81c5146f046, 2= -4313720312531179750], [1= 2, 2= 2], [1= 3, 2= 3], [1= 4, 2= 4], [1= 5, 2= 5], [1= 6, 2= 6], [1= 7, 2= 7], [1= 8, 2= 8], [1= 9, 2= 9], [1= 0, 2= 0], [1= 11, 2= 11], [1= 12, 2= 12], [1= 13, 2= 13], [1= 14, 2= 14], [1= 15, 2= 15], [1= addrof ptr $__muid_func_def_tab$$Arith_jar, 2= addrof ptr $__muid_func_def_tab$$Arith_jar], [1= 0, 2= 0], [1= addrof ptr $__muid_func_inf_tab$$Arith_jar, 2= addrof ptr $__muid_func_inf_tab$$Arith_jar], [1= addrof ptr $__muid_func_undef_tab$$Arith_jar, 2= addrof ptr $__muid_func_undef_tab$$Arith_jar], [1= addrof ptr $__muid_data_def_tab$$Arith_jar, 2= addrof ptr $__muid_data_def_tab$$Arith_jar], [1= 0, 2= 0], [1= addrof ptr $__muid_data_undef_tab$$Arith_jar, 2= addrof ptr $__muid_data_undef_tab$$Arith_jar], [1= addrof ptr $__muid_func_def_muid_tab$$Arith_jar, 2= addrof ptr $__muid_func_def_muid_tab$$Arith_jar], [1= addrof ptr $__muid_func_undef_muid_tab$$Arith_jar, 2= addrof ptr $__muid_func_undef_muid_tab$$Arith_jar], [1= addrof ptr $__muid_data_def_muid_tab$$Arith_jar, 2= addrof ptr $__muid_data_def_muid_tab$$Arith_jar], [1= addrof ptr $__muid_data_undef_muid_tab$$Arith_jar, 2= addrof ptr $__muid_data_undef_muid_tab$$Arith_jar], [1= addrof ptr $__muid_func_muid_idx_tab$$Arith_jar, 2= addrof ptr $__muid_func_muid_idx_tab$$Arith_jar], [1= 0, 2= 0]]
var $__reflection_run_hot_strtab$$Arith_jar fstatic <[26] u8> = [0, 76, 65, 114, 105, 116, 104, 59, 0, 114, 97, 105, 115, 101, 95, 115, 105, 103, 102, 112, 101, 124, 40, 41, 86, 0]
var $__reflection_strtab$$Arith_jar fstatic <[71] u8> = [0, 114, 97, 105, 115, 101, 95, 115, 105, 103, 102, 112, 101, 0, 40, 41, 86, 0, 48, 33, 48, 0, 109, 97, 105, 110, 0, 40, 91, 76, 106, 97, 118, 97, 47, 108, 97, 110, 103, 47, 83, 116, 114, 105, 110, 103, 59, 41, 86, 0, 60, 105, 110, 105, 116, 62, 0, 84, 101, 115, 116, 77, 97, 105, 110, 0, 40, 73, 41, 73, 0]
var $__compilerVersionNum$$Arith_jar <[0] <* void>> = [1, 0]
func &LArith_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LArith_3B>>) void {
  funcid 48153
  var %Reg1_R43694 <* <$LArith_3B>> localrefvar
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>> localrefvar
  var %__muid_symptr <* void>

  intrinsiccall MCCIncRef (dread ref %_this)
  intrinsiccall MCCDecRef (dread ref %Reg1_R43694)
  dassign %Reg1_R43694 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  regassign ptr %1 (dread ref %Reg1_R57)
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43694))
  intrinsiccall MCCIncRef (dread ref %Reg1_R57)
  intrinsiccall MCCDecRef (regread ptr %1)
  #Call function:Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V
  dassign %__muid_symptr 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[13] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_func_undef_tab$$Arith_jar, constval i64 9)))
  icallassigned (dread ptr %__muid_symptr, dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  intrinsiccall MPL_CLEANUP_LOCALREFVARS (dread ref %Reg1_R43694, dread ref %Reg1_R57)
  return ()
}
func &LArith_3B_7Craise__sigfpe_7C_28_29V private static native local () void {
  funcid 48154
  var %env_ptr <* void>

  callassigned &MCC_PreNativeCall (addrof ptr $__cinf_LArith_3B) { dassign %env_ptr 0 }
  regassign ptr %1 (iread ptr <* <* void>> 0 (array 0 ptr <* <[1] <* void>>> (addrof ptr $__reg_jni_func_tab$$Arith_jar, constval i32 0)))
  regassign ptr %2 (lshr ptr (regread ptr %1, constval u64 56))
  if (eq u1 ptr (regread ptr %2, constval ptr 255)) {
    callassigned &MCC_FindNativeMethodPtrWithoutException (array 0 ptr <* <[1] <* void>>> (addrof ptr $__reg_jni_func_tab$$Arith_jar, constval i32 0)) { regassign ptr %1}
    callassigned &MCC_DummyNativeMethodPtr () { regassign ptr %3}
    if (eq u1 ptr (regread ptr %1, regread ptr %3)) {
      regassign ptr %1 (addroffunc ptr &Java_Arith_raise_1sigfpe__)
      iassign <* <* void>> 0 (
        array 0 ptr <* <[1] <* void>>> (addrof ptr $__reg_jni_func_tab$$Arith_jar, constval i32 0),
        regread ptr %1)
    }
  }
  call &MCC_CallSlowNative2 (regread ptr %1, dread ptr %env_ptr, addrof ptr $__cinf_LArith_3B)
  callassigned &MCC_PostNativeCall (dread ptr %env_ptr) {}
  callassigned &MCC_CheckThrowPendingException () {}
  callassigned &MCC_SetReliableUnwindContext () {}
}
func &LArith_3B_7CTestMain_7C_28I_29I public static (var %Reg4_I i32) i32 {
  funcid 48155
  var %Reg0_I i32
  var %Reg2_I i32
  var %Reg0_R19929 <* <$Ljava_2Flang_2FArithmeticException_3B>> localrefvar
  var %Reg0_R510 <* <$Ljava_2Flang_2FThrowable_3B>> localrefvar
  var %Reg3_R19929 <* <$Ljava_2Flang_2FArithmeticException_3B>> localrefvar
  var %__muid_symptr <* void>

  #INSTIDX : 0||0000:  iconst_0
  dassign %Reg0_I 0 (constval i32 0)
  #INSTIDX : 1||0001:  istore_1
  dassign %Reg2_I 0 (dread i32 %Reg0_I)
  try { @label2 }
  #INSTIDX : 2||0002:  iinc
  dassign %Reg2_I 0 (add i32 (dread i32 %Reg2_I, constval i32 1))
  #INSTIDX : 5||0005:  iload_0
  #INSTIDX : 6||0006:  ifle
  brtrue @label0 (le i32 i32 (dread i32 %Reg4_I, constval i32 0))
  #INSTIDX : 9||0009:  invokestatic
  #Call function:LArith_3B_7Craise__sigfpe_7C_28_29V
  dassign %__muid_symptr 0 (iread ptr <* <$MUIDFuncDefTabEntry>> 1 (array 0 ptr <* <[4] <$MUIDFuncDefTabEntry>>> (addrof ptr $__muid_func_def_tab$$Arith_jar, constval i64 1)))
  icallassigned (dread ptr %__muid_symptr) {}
  #INSTIDX : 12||000c:  goto
  goto @label1
@label0   #INSTIDX : 15||000f:  iload_0
  #INSTIDX : 16||0010:  ifne
  brtrue @label1 (ne i32 i32 (dread i32 %Reg4_I, constval i32 0))
  #INSTIDX : 19||0013:  iinc
  dassign %Reg2_I 0 (add i32 (dread i32 %Reg2_I, constval i32 103))
  #INSTIDX : 22||0016:  new
  regassign ptr %1 (dread ref %Reg0_R19929)
  dassign %Reg0_R19929 0 (gcmalloc ref <$Ljava_2Flang_2FArithmeticException_3B>)
  intrinsiccall MCCDecRef (regread ptr %1)
  #INSTIDX : 25||0019:  dup
  #INSTIDX : 26||001a:  invokespecial
  #Call function:Ljava_2Flang_2FArithmeticException_3B_7C_3Cinit_3E_7C_28_29V
  dassign %__muid_symptr 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[13] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_func_undef_tab$$Arith_jar, constval i64 10)))
  icallassigned (dread ptr %__muid_symptr, dread ref %Reg0_R19929) {}
  #INSTIDX : 29||001d:  athrow
  regassign ptr %2 (dread ref %Reg0_R510)
  dassign %Reg0_R510 0 (retype ref <* <$Ljava_2Flang_2FThrowable_3B>> (dread ref %Reg0_R19929))
  intrinsiccall MCCIncRef (dread ref %Reg0_R510)
  intrinsiccall MCCDecRef (regread ptr %2)
  throw (dread ref %Reg0_R510)
@label1   #INSTIDX : 30||001e:  iinc
  dassign %Reg2_I 0 (add i32 (dread i32 %Reg2_I, constval i32 3))
  endtry
  dassign %Reg2_I 0 (cvt i32 i32 (dread i32 %Reg2_I))
  #INSTIDX : 33||0021:  goto
  goto @label3
@label2   catch { <* <$Ljava_2Flang_2FArithmeticException_3B>> }
  intrinsiccall MCCDecRef (dread ref %Reg0_R19929)
  dassign %Reg0_R19929 0 (regread ptr %%thrownval)
  #INSTIDX : 36||0024:  astore_2
  intrinsiccall MCCIncRef (dread ref %Reg0_R19929)
  intrinsiccall MCCDecRef (dread ref %Reg3_R19929)
  dassign %Reg3_R19929 0 (dread ref %Reg0_R19929)
  #INSTIDX : 37||0025:  iinc
  dassign %Reg2_I 0 (add i32 (dread i32 %Reg2_I, constval i32 100))
  dassign %Reg2_I 0 (cvt i32 i32 (dread i32 %Reg2_I))
@label3   #INSTIDX : 40||0028:  iload_1
  #INSTIDX : 41||0029:  ireturn
  intrinsiccall MPL_CLEANUP_LOCALREFVARS (dread ref %Reg0_R19929, dread ref %Reg0_R510, dread ref %Reg3_R19929)
  return (dread i32 %Reg2_I)
}
func &LArith_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg2_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void {
  funcid 48156
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>> localrefvar
  var %Reg1_I i32
  var %__muid_symptr <* void>

  #INSTIDX : 0||0000:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %1 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[4] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$Arith_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %1)
  #INSTIDX : 3||0003:  bipush
  dassign %Reg1_I 0 (constval i32 -5)
  #INSTIDX : 5||0005:  invokestatic
  #Call function:LArith_3B_7CTestMain_7C_28I_29I
  dassign %__muid_symptr 0 (iread ptr <* <$MUIDFuncDefTabEntry>> 1 (array 0 ptr <* <[4] <$MUIDFuncDefTabEntry>>> (addrof ptr $__muid_func_def_tab$$Arith_jar, constval i64 2)))
  icallassigned (dread ptr %__muid_symptr, dread i32 %Reg1_I) { dassign %Reg1_I 0 }
  #INSTIDX : 8||0008:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 272)),
    dread ref %Reg0_R562,
    dread i32 %Reg1_I) {}
  #INSTIDX : 11||000b:  getstatic
  intrinsiccall MPL_CLINIT_CHECK (addrof ptr $__cinf_Ljava_2Flang_2FSystem_3B)
  regassign ptr %2 (dread ref %Reg0_R562)
  #Read from: Ljava_2Flang_2FSystem_3B_7Cout
  dassign %Reg0_R562 0 (iread ref <* <* <$Ljava_2Fio_2FPrintStream_3B>>> 0 (iread ptr <* <$MUIDUnifiedUndefTabEntry>> 1 (array 0 ptr <* <[4] <$MUIDUnifiedUndefTabEntry>>> (addrof ptr $__muid_data_undef_tab$$Arith_jar, constval i64 1))))
  intrinsiccall MCCIncRef (dread ref %Reg0_R562)
  intrinsiccall MCCDecRef (regread ptr %2)
  #INSTIDX : 14||000e:  iconst_0
  dassign %Reg1_I 0 (constval i32 0)
  #INSTIDX : 15||000f:  invokestatic
  #Call function:LArith_3B_7CTestMain_7C_28I_29I
  dassign %__muid_symptr 0 (iread ptr <* <$MUIDFuncDefTabEntry>> 1 (array 0 ptr <* <[4] <$MUIDFuncDefTabEntry>>> (addrof ptr $__muid_func_def_tab$$Arith_jar, constval i64 2)))
  icallassigned (dread ptr %__muid_symptr, dread i32 %Reg1_I) { dassign %Reg1_I 0 }
  #INSTIDX : 18||0012:  invokevirtual
  icallassigned (
    iread u64 <* u64> 0 (add ptr (
      iread ptr <* <$__class_meta__>> 6 (iread ref <* <$Ljava_2Flang_2FObject_3B>> 1 (dread ref %Reg0_R562)),
      constval u32 272)),
    dread ref %Reg0_R562,
    dread i32 %Reg1_I) {}
  #INSTIDX : 21||0015:  return
  intrinsiccall MPL_CLEANUP_LOCALREFVARS (dread ref %Reg0_R562)
  return ()
}
func &Java_Arith_raise_1sigfpe__ weak () void {
  funcid 0

  callassigned &MCC_CannotFindNativeMethod (conststr ptr "LArith_3B_7Craise__sigfpe_7C_28_29V") {}
}

max_stars_repo_name: Bravewtz/OpenArkCompiler

max_issues_repo_name: Bravewtz/OpenArkCompiler

max_forks_repo_name: Bravewtz/OpenArkCompiler


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: lda_exc *)
(* prefix:
  lda_c_pw_params *params;

  assert(p->params != NULL);
  params = (lda_c_pw_params * )(p->params);
*)

$ifdef lda_c_pw_params
params_a_pp     := [1,  1,  1]:
params_a_a      := [0.031091, 0.015545, 0.016887]:
params_a_alpha1 := [0.21370,  0.20548,  0.11125]:
params_a_beta1  := [7.5957, 14.1189, 10.357]:
params_a_beta2  := [3.5876, 6.1977, 3.6231]:
params_a_beta3  := [1.6382, 3.3662,  0.88026]:
params_a_beta4  := [0.49294, 0.62517, 0.49671]:
params_a_fz20   := 1.709921:
$endif

$ifdef lda_c_pw_modified_params
params_a_a      := [0.0310907, 0.01554535, 0.0168869]:
params_a_fz20   := 1.709920934161365617563962776245:
$endif

(* Equation (10) *)
g_aux := (k, rs) -> params_a_beta1[k]*sqrt(rs) + params_a_beta2[k]*rs
  + params_a_beta3[k]*rs^1.5 + params_a_beta4[k]*rs^(params_a_pp[k] + 1):
g     := (k, rs) -> -2*params_a_a[k]*(1 + params_a_alpha1[k]*rs)
  * log(1 +  1/(2*params_a_a[k]*g_aux(k, rs))):

(* Equation (8) *)
(* Attention, the function g parametrizes -alpha *)
f_pw := (rs, zeta) ->
  g(1, rs) + zeta^4*f_zeta(zeta)*(g(2, rs) - g(1, rs) + g(3, rs)/params_a_fz20)
  - f_zeta(zeta)*g(3, rs)/params_a_fz20:

f := (rs, zeta) -> f_pw(rs, zeta):


text:
############################################
#  Performing finite differencing  using
#  FD toolkit. See previous script
#  'compactification.mpl' for how the tt table
#  is generated
############################################

read "/home/arman/FD/FD.mpl":
CFD();
MFD();

interface(warnlevel=0);
with(StringTools):

##############################################################
# By default FD uses second order centered in Gen_Sten for RHS
# which is exactly what we want
##############################################################

for ii from 1 to nops(all_tensors) do
  nm:= all_tensors[ii][1]:
  rank := all_tensors[ii][2]:
  tp  := all_tensors[ii][3]:

  printf("Discretizing %a\n",nm);
  if rank <> 0 then
      # Note the collect operator on powers of 1/hx, it will help the compiler's O3 flag a lot!
      tt[nm][dis] := table( [ seq(  coord[jj] = collect(Gen_Sten(tt[nm][all][coord[jj]]),1/hx) ,jj=1..2 ) ] );
      tt[nm][disreg] := table( [seq( coord[jj] = collect(Gen_Sten(tt[nm][reg][coord[jj]]),1/hx) ,jj=1..2 ) ] );
  else
      tt[nm][dis] := collect(Gen_Sten(tt[nm][all]),1/hx);
      tt[nm][disreg] := collect(Gen_Sten(tt[nm][reg]),1/hx);
  end if;

end do:

###################################################
# Bounday condition at infinity for the tensors:
# Note that we use 'myzero*x' to indicate '0' grid
# function at infinity, See documentation of FD
# for why this is neccessary
# The overhead of this is almost zero!
###################################################

inf_b_cond:= table( [
EMPH = myzero*x+1,
DPH = myzero*x,
DB  = myzero*x,
dg  = [myzero*x,myzero*x],
met = [myzero*x+1,myzero*x+1],
DDL = [myzero*x,myzero*x],
DTK = myzero*x,
RR  = [myzero*x,myzero*x],
U   = [myzero*x,myzero*x],
P   = [myzero*x,myzero*x],
DAA = [myzero*x,myzero*x],
DLAM = [myzero*x,myzero*x],
C1 = myzero*x,
C5 = myzero*x,
HS = psi(t,x)-1 + myzero*x,
HS3 = psi(t,x)-1 + myzero*x,
EQRPIb2 = myzero*x,
EQIPIb2 = myzero*x,
EQIPSI = myzero*x,
EQRPSI = myzero*x,
DRPSI = myzero*x,
DIPSI = myzero*x
]);

# Discretizing the infinity boundary conditions

for ii from 1 to nops(all_tensors) do
  nm:= all_tensors[ii][1]:
  rank := all_tensors[ii][2]:
  tp  := all_tensors[ii][3]:

  printf("Discretizing %a\n",nm);
  if rank <> 0 then
      tt[nm][disinf] := table( [ seq(  coord[jj] = Gen_Sten(inf_b_cond[nm][jj]) ,jj=1..2 ) ] );
  else
      tt[nm][disinf] := Gen_Sten(inf_b_cond[nm]);
  end if;

end do:

# Updating FD's table to default centered in all coordinates

pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

############################################
# A table to hold the LHS of BSSN equations
# in their continuous form
############################################

evol_eq_LHS_tbl := table ( [
DPH = diff(phi(t,x),t),
dg = [diff(A(t,x),t),diff(B(t,x),t)],
DTK = diff(K(t,x),t),
DAA = [diff(Axx(t,x),t),diff(Athth(t,x),t)],
DLAM = [diff(Lamx(t,x),t),0],
EQRPIb2 = diff(rPIb2(t,x),t),
EQIPIb2 = diff(iPIb2(t,x),t),
EQRPSI  = diff(rpsi(t,x),t),
EQIPSI  = diff(ipsi(t,x),t)
] );

###############################################
# A table to hold what each tensorial equation
# is actually will be solved for, this will be
# needed by FD to decide what to solve the FDA
# for!
###############################################

solve_for_tbl := table ( [
   DPH = phi(n+1,i),
   dg  = table( [ x = A(n+1,i) , theta=B(n+1,i)  ] ),
   DTK = K(n+1,i),
   DAA = table( [ x = Axx(n+1,i), theta=Athth(n+1,i) ] ),
   DLAM = table( [ x = Lamx(n+1,i), theta = 0 ] ),
   EQRPIb2 = rPIb2(n+1,i),
   EQIPIb2 = iPIb2(n+1,i),
   EQRPSI = rpsi(n+1,i),
   EQIPSI = ipsi(n+1,i)
] );

#################################################
# A table to hold what are the actual discrete
# form of the functions are, this will be needed
# by FD
#################################################

dis_func_tbl := table ( [
   DPH = phi(n,i),
   dg  = table( [ x = A(n,i) , theta=B(n,i)  ] ),
   DTK = K(n,i),
   DAA = table( [ x = Axx(n,i), theta=Athth(n,i) ] ),
   DLAM = table( [ x = Lamx(n,i), theta = 0 ] ),
   EQRPIb2 = rPIb2(n,i),
   EQIPIb2 = iPIb2(n,i),
   EQRPSI = rpsi(n,i),
   EQIPSI = ipsi(n,i)
] );

##########################################################
# The implementaion of the famous Kreiss/Oliger disspation
# using FD's 'FD' operator, will be used in the RHS of the
# evolution equations to suppress numerical noise
# See FD's documentation on how 'FD' operator works, but
# if you are familiar with RNPL, it should be pretty clear
# what is happening here!
##########################################################

DISS := f-> -zepsdis / (16*ht) * ( 6*FD(f,[[0],[0,0,0]])
                                   + FD(f,[[0],[2,0,0]])
                                   + FD(f,[[0],[-2,0,0]])
                                 - 4*FD(f,[[0],[1,0,0]])
                                 - 4*FD(f,[[0],[-1,0,0]])
                                );

# Defining average operator for Crank-Nickelson second order update
AVGT := f -> (FD(f,[[1],[0,0,0]]) + FD(f,[[0],[0,0,0]]))/2 ;

# Advanced time operator
NT := f -> FD(f,[[1],[0,0,0]]);

# 2nd order interpolation at the boundary (used at origin)
QFIT := f-> FD(f,[[0],[0,0,0]]) - 4*FD(f,[[0],[1,0,0]])/3 + FD(f,[[0],[2,0,0]])/3 ;

######################################################################
# Shift operator, will be used for the point next to origin condition
# See Eq. (58) in the paper: http://arxiv.org/pdf/1508.01614v2.pdf
#######################################################################

SB := f -> FD(f,[[0],[-1,0,0]]);

# Updating FD's table to default in all coordinates
pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

# The evolutionary type equations:
evol_trm := [[DPH,0],[dg,2],[DTK,0],[DAA,2],[DLAM,1], [EQRPIb2,0] , [EQIPIb2,0], [EQRPSI,0], [EQIPSI,0] ];

# Initiating a table to hold the actual equations (residuals) in discrete form
evol_eq_tbl := table ( [] );

# Updating to forward in time differencing scheme, i.e. Crank-Niecklson when combined with AVGT
FD_table[t] := [ [0] , [0,1] ];

# List of all even/odd functions, will be used by FD for imposing inner-boundary conditions
set_of_even_funcs:={alpha,K,phi,Axx,Athth,A,B,psi,a1,a2,b1,ipsi,rpsi,iPIb2,rPIb2,rPI,iPI,divbeta}:
set_of_odd_funcs:={beta,Lamx,rPHI,iPHI}:

for ii from 1 to nops(evol_trm) do
  nm := evol_trm[ii][1];
  rank := evol_trm[ii][2];
  if rank <> 0 then
    for jj from 1 to 2 do
       cn := coord[jj];
       #######  THIS IS THE EQUATION IN FDA FORM ############
       evol_eq_tbl[nm][dis][cn] := Gen_Sten(evol_eq_LHS_tbl[nm][jj]) - AVGT(tt[nm][dis][cn]);
       evol_eq_tbl[nm][disreg][cn] := Gen_Sten(evol_eq_LHS_tbl[nm][jj]) - AVGT(tt[nm][disreg][cn]);
       # Applying boundary conditions using FD's facilities at the origin
       t_mp_1:= A_FD_Even(tt[nm][disreg][cn],x,set_of_even_funcs,0,"forward");
       t_mp_2:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");
       # Storing  (in 'disorg') a version of FDA that can be used next to the origin
       evol_eq_tbl[nm][disorg][cn] := Gen_Sten(evol_eq_LHS_tbl[nm][jj]) - AVGT(t_mp_2);
       # Storing the form of the equation at infinity:
       evol_eq_tbl[nm][disinf][cn] := Gen_Sten(evol_eq_LHS_tbl[nm][jj]) - AVGT(tt[nm][disinf][cn]);
    end do
  else
      # Same as above
      evol_eq_tbl[nm][dis] := Gen_Sten(evol_eq_LHS_tbl[nm]) - AVGT(tt[nm][dis]);
      evol_eq_tbl[nm][disreg] := Gen_Sten(evol_eq_LHS_tbl[nm]) - AVGT(tt[nm][disreg]);
      t_mp_1:= A_FD_Even(tt[nm][disreg],x,set_of_even_funcs,0,"forward");
      t_mp_2:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");
      evol_eq_tbl[nm][disorg] := Gen_Sten(evol_eq_LHS_tbl[nm]) - AVGT(t_mp_2);
      evol_eq_tbl[nm][disinf] := Gen_Sten(evol_eq_LHS_tbl[nm]) - AVGT(tt[nm][disinf]);
  end if;
end do;

# The evaluation type equations, i.e. 'work variables'
eval_term := [[EMPH,0] , [DB,0] ,[met,2],[DDL,2] ,[RR,2],[U,2], [P,2] ,[C1,0], [C5,0], [DIPSI,0], [DRPSI,0]];

eval_expr_tbl := table ( [] );

for ii from 1 to nops(eval_term) do
  nm := eval_term[ii][1];
  rank := eval_term[ii][2];
  if rank <> 0 then
    for jj from 1 to 2 do
      cn := coord[jj];
      eval_expr_tbl[nm][dis][cn] := (tt[nm][dis][cn]);
      eval_expr_tbl[nm][disreg][cn] := (tt[nm][disreg][cn]);
      eval_expr_tbl[nm][disinf][cn] := (tt[nm][disinf][cn]);
      # Same as evolution case
      t_mp_1:= A_FD_Even(tt[nm][disreg][cn],x,set_of_even_funcs,0,"forward");
      t_mp_2:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");
      eval_expr_tbl[nm][disorg][cn] := t_mp_2;
    end do;
  else
      eval_expr_tbl[nm][dis]:= (tt[nm][dis]);
      eval_expr_tbl[nm][disreg]:= (tt[nm][disreg]);
      eval_expr_tbl[nm][disinf]:= (tt[nm][disinf]);
      t_mp_1:= A_FD_Even(tt[nm][disreg],x,set_of_even_funcs,0,"forward");
      t_mp_2:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");
      eval_expr_tbl[nm][disorg] := t_mp_2;
  end if;
end do;

# Switching to forward in 'x' (ie 'r') for doing FDA at the origin
pl:=table([ t=[-1,-1],x=[0,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

# Again in time, we use CN scheme
FD_table[t] := [ [0] , [0,1] ];

#######################################################
# Table for boundary conditions at the origin
# See the paper: http://arxiv.org/pdf/1508.01614v2.pdf
#######################################################

org_bdc_tbl := table ( [
   # d phi / dx must be zero at origin
   DPH = NT(Gen_Sten(diff(phi(t,x),x))),
   # conformal metric components are 1 at the origin
   dg  = table( [ x = A(n+1,i) - 1 +myzero*x(i) , theta= B(n+1,i) -1 + myzero*x(i)  ] ),
   # exterinsic curvature has zero spatial derivative at origin:
   DTK = NT(Gen_Sten(diff(K(t,x),x))),
   # (from dg eq above) conformal extrinsic curvature are fixed 0 at the origin
   DAA = table( [ x = Axx(n+1,i) +myzero*x(i), theta= Athth(n+1,i)+myzero*x(i) ] ),
   # Conformal connection is vector thus zero at the origin
   DLAM = table( [ x = Lamx(n+1,i) +myzero*x(i), theta= 0 ] ),
   # We use interpolation for matter fields
   EQRPIb2 = NT(QFIT(rPIb2(n,i))),
   EQIPIb2 = NT(QFIT(iPIb2(n,i))),
   EQRPSI  = NT(QFIT(rpsi(n,i))),
   EQIPSI  = NT(QFIT(ipsi(n,i)))
] );

##################################################################
# Extra boundary conditions for the point next to the origin
# as explained in the paper: http://arxiv.org/pdf/1508.01614v2.pdf
# See eq. 58 in the paper
###################################################################

orgp1_bdc_tbl := table ( [
   # Just using the original version, (no extra condition)
   DPH = evol_eq_tbl[DPH][dis],
   # Using higher order derivative condition (see the paper)
   dg  = table( [ x = NT(SB(Gen_Sten(diff(A(t,x),x)))) ,theta= NT(SB(Gen_Sten(diff(B(t,x),x)))) ] ),
   # No BC here
   DTK = evol_eq_tbl[DTK][dis],
   # Same as dg
   DAA  = table( [ x = NT(SB(Gen_Sten(diff(Axx(t,x),x)))) +hx*myzero ,theta= hx*myzero+NT(SB(Gen_Sten(diff(Athth(t,x),x))))  ] ),
   # Same as dg
   DLAM = table( [ x = NT(SB(Gen_Sten(diff(Lamx(t,x),x,x))) ), theta= 0 ] ),
   # No BC here
   EQRPIb2 = evol_eq_tbl[EQRPIb2][dis],
   EQIPIb2 = evol_eq_tbl[EQIPIb2][dis],
   EQRPSI  = evol_eq_tbl[EQRPSI][dis],
   EQIPSI  = evol_eq_tbl[EQIPSI][dis]

] );

# Adding the newly defined FDA's to the original table:
for ii from 1 to nops(evol_trm) do
  nm := evol_trm[ii][1];
  rank := evol_trm[ii][2];
  if rank <> 0 then
    for jj from 1 to 2 do
       cn := coord[jj];
        evol_eq_tbl[nm][disorg][cn] := org_bdc_tbl[nm][cn];
        evol_eq_tbl[nm][disorgp1][cn] := orgp1_bdc_tbl[nm][cn];
    end do
  else
       evol_eq_tbl[nm][disorg] := org_bdc_tbl[nm];
       evol_eq_tbl[nm][disorgp1] := orgp1_bdc_tbl[nm];
  end if;
end do:

# BC for work variables, in reality the values will never be used
e_org_bdc_tbl := table ( [
EMPH = myzero*x+1,
met = [myzero*x+1,myzero*x+1],
DDL = [myzero*x,myzero*x],
RR  = [myzero*x,myzero*x],
U   = [myzero*x,myzero*x],
P   = [myzero*x,myzero*x],
C1 = myzero*x,
C5 = myzero*x,
DRPSI = myzero*x,
DIPSI = myzero*x,
DB  = myzero*x
 ] );

pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

# Adding them to the equation table:
for ii from 1 to nops(eval_term) do
  nm := eval_term[ii][1];
  rank := eval_term[ii][2];
  if rank <> 0 then
    for jj from 1 to 2 do
       cn := coord[jj];
        eval_expr_tbl[nm][disorg][cn] := Gen_Sten(e_org_bdc_tbl[nm][jj]);
    end do
  else
        eval_expr_tbl[nm][disorg] := Gen_Sten(e_org_bdc_tbl[nm]);
  end if;
end do:

# Treating Hamiltonian constraint seperately
pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

nm := HS;

const_eq_tbl := table ( [] );

const_eq_tbl[nm][dis] := tt[nm][dis];
const_eq_tbl[nm][disreg] := tt[nm][disreg];
t_mp_1 := A_FD_Even(tt[nm][disreg],x,set_of_even_funcs,0,"forward");
t_mp_2 := A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");
const_eq_tbl[nm][disorg] := t_mp_2;

pl:=table([ t=[-1,-1],x=[-1,0],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

# BC for Hamiltonian constraint at infinity:
const_eq_tbl[nm][inf] := Gen_Sten(ctfm(x)*diff(psi(t,x),x)/ctfmp(x) + psi(t,x) - 1);

pl:=table([ t=[-1,-1],x=[0,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

const_eq_tbl[nm][disorg] := Gen_Sten(diff(psi(t,x),x));

pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

# Will store header files, to pre-generate the C routine
Other_Headers := NULL;

# DDS for HC is easy:

spec := [
  { i = [1,1,1],  b=xmin } =  const_eq_tbl[nm][disorg] ,
  { i = [2,Nx-1,1] } = const_eq_tbl[nm][dis] ,
  { i = [Nx,Nx,1], b=xmax  } = const_eq_tbl[nm][inf]
 ];

# Residual evaluator for HC
A_Gen_Res_Code(spec,input="d",proc_name="resid_hs");
Other_Headers := Other_Headers , "resid_hs";

# Solver for HC
A_Gen_Solve_Code(spec,{psi(n,i)},proc_name="solve_hs");
Other_Headers := Other_Headers , "solve_hs";

# Now working with momentum constraint
pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

nm := MK;

const_eq_tbl := table ( [] );

const_eq_tbl[nm][dis] := tt[nm][dis][x];
const_eq_tbl[nm][disreg] := tt[nm][disreg][x];
t_mp_1 := A_FD_Even(tt[nm][disreg][x],x,set_of_even_funcs,0,"forward");
t_mp_2 := A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");
const_eq_tbl[nm][disorg] := t_mp_2;

spec := [
  { i = [1,1,1],  b=xmin } =  const_eq_tbl[nm][disorg] ,
  { i = [2,Nx-1,1] } = const_eq_tbl[nm][dis] ,
  { i = [Nx,Nx,1], b=xmax  } = myzero*x(i)
 ];

# Evaluator routine for momentum constraint, for testing purposes
A_Gen_Eval_Code(spec,input="d",proc_name="ire_mk");

##########################################################
# See the bssn_spher.mw for what HS3 is, another form of
# Hamiltonian constraint
##########################################################
pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

nm := HS3;

const_eq_tbl := table ( [] );

const_eq_tbl[nm][dis] := tt[nm][dis];

pl:=table([ t=[-1,-1],x=[-1,0],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl);

const_eq_tbl[nm][inf] := Gen_Sten(psi(t,x) - 1 +myzero*ctfm(x));

pl:=table([ t=[-1,-1],x=[0,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

const_eq_tbl[nm][disorg] := Gen_Sten(diff(psi(t,x),x)):

pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

Other_Headers := NULL;

spec := [
  { i = [1,1,1],  b=xmin } =  const_eq_tbl[nm][disorg] ,
  { i = [2,Nx-1,1] } = const_eq_tbl[nm][dis] ,
  { i = [Nx,Nx,1], b=xmax  } = const_eq_tbl[nm][inf]
 ]:

A_Gen_Res_Code(spec,input="d",proc_name="resid_hs3");
Other_Headers := Other_Headers , "resid_HS3";

A_Gen_Solve_Code(spec,{psi(n,i)},proc_name="solve_hs3");
Other_Headers := Other_Headers , "solve_HS3";

# End of constraint equation solvers / evaluators

# Initializer expressions for grid functions:

init_alpha_expr:=myzero*x+1;
init_beta_expr:=myzero*x;
init_A_expr:=myzero*x+1;
init_B_expr:=myzero*x+1;
init_Axx_expr:=myzero*x;
init_Athth_expr:=myzero*x;
init_TK_expr:=myzero*x;
init_Lamx_expr:=myzero*x;

# Some redefintion of variables:

phi_expr:=log(psi(t,x));
psi_expr := exp(phi(t,x));
a2_expr:= psi(t,x)^4*A(t,x);
a1_expr := sqrt(a2(t,x));
b1_expr := psi(t,x)^2*sqrt(B(t,x));

lv := 'lv';
rv := 'rv';

# Initializer for matter field (Gaussian pulse on real coordinate)
rpsi_expr :=  amp*exp( - ( ctfm(x) - x0 )^2 / deltx^2 );
ipsi_expr := myzero*x*x;

# Note the correction due to compactification
init_rPHI_expr := diff(rpsi_expr, x)/ctfmp(x);
init_iPHI_expr := diff(ipsi_expr, x)/ctfmp(x);

rpsidot_expr := myzero*x;
ipsidot_expr := omega*rpsi_expr;

# PI = a/alpha(psidot - beta*psiprime)

init_rPI_expr := a1(t,x)/alpha(t,x)*(rpsidot(t,x) - beta(t,x)*rPHI(t,x));
init_iPI_expr := a1(t,x)/alpha(t,x)*(ipsidot(t,x) - beta(t,x)*iPHI(t,x));

init_rPIb2_expr := rPI(t,x)*b1(t,x)^2;
init_iPIb2_expr := iPI(t,x)*b1(t,x)^2;

rPI_comp_expr :=  rPIb2(t,x)/b1(t,x)^2;
iPI_comp_expr :=  iPIb2(t,x)/b1(t,x)^2;

PI2_expr := (rPI(t,x)^2 + iPI(t,x)^2);
PHI2_expr := (rPHI(t,x)^2 + iPHI(t,x)^2);

# Energy-momentum tensor expressions:

U_expr := mass*(rpsi(t,x)^2 + ipsi(t,x)^2);

rho_expr:= ( PI2(t,x) + PHI2(t,x) ) / ( 2*a2(t,x)) + U(t,x)/2;

JSxx_expr:=  ( PI2(t,x) + PHI2(t,x) ) / (2) - a2(t,x)*U(t,x)/2;

JSthetatheta_expr :=  (  (PI2(t,x) - PHI2(t,x)) / ( 2*a2(t,x) ) - U(t,x)/2  ) * b1(t,x)^2;

TS_expr := ( 3*PI2(t,x) - PHI2(t,x) ) / (2*a2(t,x) ) - 3/2*U(t,x) ;
Sx_expr := - (  rPI(t,x)*rPHI(t,x) + iPI(t,x)*iPHI(t,x)  ) / (a1(t,x))  ;

############################################################
# all_expr will be used for pre-generation of main.c routine
# this is an experimental feature in FD
############################################################
all_expr := rPHI_expr+mass+omega;

# Creating initializer routines

Gen_Eval_Code(init_alpha_expr,input="c",proc_name="init_alpha");
Gen_Eval_Code(init_beta_expr,input="c",proc_name="init_beta");

Gen_Eval_Code(init_A_expr,input="c",proc_name="init_a");
Gen_Eval_Code(init_B_expr,input="c",proc_name="init_b");

Gen_Eval_Code(init_TK_expr,input="c",proc_name="init_tk");
Gen_Eval_Code(init_Axx_expr,input="c",proc_name="init_axx");
Gen_Eval_Code(init_Athth_expr,input="c",proc_name="init_athth");
Gen_Eval_Code(init_Lamx_expr,input="c",proc_name="init_lamx");

Gen_Eval_Code(rpsi_expr,input="c",proc_name="init_rpsi");
Gen_Eval_Code(ipsi_expr,input="c",proc_name="init_ipsi");

Gen_Eval_Code(init_iPHI_expr,input="c",proc_name="init_iphi");
Gen_Eval_Code(init_rPHI_expr,input="c",proc_name="init_rphi");

Gen_Eval_Code(rpsidot_expr,input="c",proc_name="init_rpsidot");
Gen_Eval_Code(ipsidot_expr,input="c",proc_name="init_ipsidot");

Gen_Eval_Code(init_rPI_expr,input="c",proc_name="init_rpi");
Gen_Eval_Code(init_iPI_expr,input="c",proc_name="init_ipi");

Gen_Eval_Code(init_iPIb2_expr,input="c",proc_name="init_ipib2");
Gen_Eval_Code(init_rPIb2_expr,input="c",proc_name="init_rpib2");

Gen_Eval_Code(phi_expr,input="c",proc_name="compute_phi");
Gen_Eval_Code(psi_expr,input="c",proc_name="compute_psi");
Gen_Eval_Code(a2_expr,input="c",proc_name="compute_a2");
Gen_Eval_Code(a1_expr,input="c",proc_name="compute_a1");
Gen_Eval_Code(b1_expr,input="c",proc_name="compute_b1");

Gen_Eval_Code(iPI_comp_expr,input="c",proc_name="compute_ipi");
Gen_Eval_Code(rPI_comp_expr,input="c",proc_name="compute_rpi");

Gen_Eval_Code(PHI2_expr,input="c",proc_name="compute_phi2");
Gen_Eval_Code(PI2_expr,input="c",proc_name="compute_pi2");

Gen_Eval_Code(U_expr,input="c",proc_name="compute_u");

Gen_Eval_Code(rho_expr,input="c",proc_name="compute_rho");
Gen_Eval_Code(JSxx_expr,input="c",proc_name="compute_jsxx");
Gen_Eval_Code(JSthetatheta_expr,input="c",proc_name="compute_jsthth");
Gen_Eval_Code(TS_expr,input="c",proc_name="compute_ts");
Gen_Eval_Code(Sx_expr,input="c",proc_name="compute_sx");

C5psi_expr := C5s(t,x)*psi(t,x)^5;
Gen_Eval_Code(C5psi_expr,input="c",proc_name="compute_c5psi");

# Newly generated header files, will be passed to C code generator
Other_Headers := Other_Headers , "init_alpha", "init_beta", "init_A" , "init_B", "init_TK", "init_Axx", "init_Athth",
                "init_Lamx", "compute_phi", "compute_psi", "compute_a2", "compute_a1",
                "compute_b1", "init_rpsi", "init_ipsi", "init_iphi", "init_rphi",
                "init_ipib2", "init_rpib2", "compute_ipi", "compute_rpi",
                "compute_phi2", "compute_pi2", "compute_u", "compute_rho", "compute_jsxx",
                "compute_jsthth", "compute_ts", "compute_sx", "init_rpsidot", "init_ipsidot", "init_rpi","init_ipi";

# Creating DDS for work variables and calculator routines
all_work_tensors:=[
      [EMPH,0,dn], [DB,0,dn],
         [met,2,dn], [DDL,2,dn],
         [RR,2,dn], [U,2,dn] , [P,2,dn],
         [C1,0,dn], [C5,0,dn], [DRPSI,0,dn], [DIPSI,0,dn]
             ];

for ii from 1 to nops(all_work_tensors) do
  nm:= all_work_tensors[ii][1]:
  rank := all_work_tensors[ii][2]:
  if rank <> 0 then
     for jj from 1 to 2 do
       cn := coord[jj];
       # DDS for work variables
       spec := [
                     { i = [1,1,1] , b=xmin }     = eval_expr_tbl[nm][disorg][cn]    ,
                     { i = [2,Nx-1,1] }  = eval_expr_tbl[nm][dis][cn]  ,
                     { i = [Nx,Nx,1] , b=xmax }  = eval_expr_tbl[nm][disinf][cn]
               ];
      # creating calculator routines fortran routines for work variables
      A_Gen_Eval_Code(spec,input="d",proc_name=LowerCase(cat("compute_",nm,cn)));
      Other_Headers := Other_Headers , LowerCase(cat("compute_",nm,cn));
     end do;
  else
        spec := [
                     { i = [1,1,1] , b=xmin }     = eval_expr_tbl[nm][disorg]    ,
                     { i = [2,Nx-1,1] }  = eval_expr_tbl[nm][dis]  ,
                     { i = [Nx,Nx,1] , b=xmax  }  = eval_expr_tbl[nm][disinf]
               ];
      A_Gen_Eval_Code(spec,input="d",proc_name=LowerCase(cat("compute_",nm)));
      Other_Headers := Other_Headers , LowerCase(cat("compute_",nm));
  end if;
end do:

# Treating conformal connection equation first:

spec := [
  { i = [1,1,1] , b=xmin } = myzero*x(i),
  { i = [2,Nx-1,1] } = tt[DLAM][dis][x] ,
  { i = [Nx,Nx,1] , b=xmax  } = myzero*x(i)
];

A_Gen_Eval_Code(spec, input="d",proc_name="eval_dlam");

# Creating a table to hold the K/O dissipation at various points of discrete domain
dis_tbl := table([]);
for ii from 1 to nops(evol_trm) do
  nm := evol_trm[ii][1];
  rank := evol_trm[ii][2];
  if rank <> 0 then
    for jj from 1 to 2 do
       cn := coord[jj];
       fn := dis_func_tbl[nm][cn];

       dis_tbl[nm][mid][cn] := -AVGT(DISS(fn));
       dis_tbl[nm][org][cn] := -AVGT(DISS(fn));
       dis_tbl[nm][orgp1][cn] :=  -AVGT(DISS(fn));

       # Using FD's even odd facility to built the correct version of FDA
       # That can be used at the origin / infinity
       t_mp_1:= A_FD_Even(dis_tbl[nm][org][cn],x,set_of_even_funcs,0,"forward");
       dis_tbl[nm][org][cn]:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");

       t_mp_1:= A_FD_Even(dis_tbl[nm][orgp1][cn],x,set_of_even_funcs,-1,"forward");
       dis_tbl[nm][orgp1][cn]:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,-1,"forward");

    end do
  else
      fn := dis_func_tbl[nm];

       dis_tbl[nm][mid] := -AVGT(DISS(fn));
       dis_tbl[nm][org] := -AVGT(DISS(fn));
       dis_tbl[nm][orgp1] :=  -AVGT(DISS(fn));

       t_mp_1:= A_FD_Even(dis_tbl[nm][org],x,set_of_even_funcs,0,"forward");
       dis_tbl[nm][org]:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,0,"forward");

       t_mp_1:= A_FD_Even(dis_tbl[nm][orgp1],x,set_of_even_funcs,-1,"forward");
       dis_tbl[nm][orgp1]:= A_FD_Odd(t_mp_1,x,set_of_odd_funcs,-1,"forward");

  end if;
end do:

####### NOW EVOLUTION EQUATIONS ##########

for ii from 1 to nops(evol_trm) do
  nm := evol_trm[ii][1]:
  rank := evol_trm[ii][2]:
  if rank <> 0 then
    for jj from 1 to 2 do
      cn := coord[jj];
      fn := dis_func_tbl[nm][cn]:
      # there are 5 conditions, origin, next to origin, middle, next to infinity and infinity
      # here is the glory of FD! (Also we are adding the right dissipation expression)
      spec := [
                 {i = [1,1,1] , b=xmin }    =  evol_eq_tbl[nm][disorg][cn],
                 {i = [2,2,1] , b=xmin }    =  evol_eq_tbl[nm][disorgp1][cn] + dis_tbl[nm][orgp1][cn],
                 {i = [3,Nx-2,1] } =  evol_eq_tbl[nm][dis][cn] + dis_tbl[nm][mid][cn],
                 {i = [Nx-1,Nx-1,1] } =  evol_eq_tbl[nm][dis][cn],
                 {i = [Nx,Nx,1] , b = xmax } =  evol_eq_tbl[nm][disinf][cn]
              ];
      if not( nm=DLAM and cn=theta) then
         print(cat("Solve code for:",nm,cn," Generated."));
         A_Gen_Solve_Code(spec,{solve_for_tbl[nm][cn] },input="d",proc_name=LowerCase(cat("evolve_",nm,cn)));
         Other_Headers := Other_Headers , cat("evolve_",nm,cn);
         A_Gen_Res_Code(spec,input="d",proc_name=LowerCase(cat("resid_",nm,cn)));
         Other_Headers := Other_Headers , cat("resid_",nm,cn);
         print(cat("Res code for:",nm,cn," Generated."));
      end if;
    end do;
  else
      fn := dis_func_tbl[nm];
      spec := [
                 {i = [1,1,1] , b=xmin }    =  evol_eq_tbl[nm][disorg],
                 {i = [2,2,1] , b=xmin }    =  evol_eq_tbl[nm][disorgp1] + dis_tbl[nm][orgp1],
                 {i = [3,Nx-2,1] } =  evol_eq_tbl[nm][dis] + dis_tbl[nm][mid],
                 {i = [Nx-1,Nx-1,1] } =  evol_eq_tbl[nm][dis],
                 {i = [Nx,Nx,1] , b=xmax  } =  evol_eq_tbl[nm][disinf]

              ];
    A_Gen_Solve_Code(spec,{solve_for_tbl[nm]},proc_name=LowerCase(cat("evolve_",nm)));
    Other_Headers := Other_Headers , cat("evolve_",nm);
    print(cat("Solve code for:",nm," Generated."));
    A_Gen_Res_Code(spec,input="d",proc_name=LowerCase(cat("resid_",nm)));
    Other_Headers := Other_Headers , cat("resid_",nm);
    print(cat("Res code for:",nm," Generated."));
  end if;
end do:

#Generating driver header:
pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):

#############################################################
# This is where I use the Other_Headers in FD
# GDC is an experimental feature, it generates a tempalte
# for a main C driver routine.
# It uses all_expr variable to learn all necessary definition
# for a solver routine, enable it and see what it generates
#############################################################

#GDC(all_expr+myzero,input="c",output="main",other_headers={Other_Headers});

# Creating IRE for BSSN equations (notice how quick this is!)

printf("Generating IREs...\n");
pl:=table([ t=[-1,-1],x=[-1,-1],y=[-1,-1],z=[-1,-1] ]):
Update_FD_Table(2,pl):
interface(warnlevel = 0):

for ii from 1 to nops(evol_trm) do
  nm := evol_trm[ii][1];
  rank := evol_trm[ii][2];
  if rank <> 0 then
    for jj from 1 to 2 do
       cn := coord[jj];
       all_expr:=all_expr + ( evol_eq_LHS_tbl[nm][jj] -  tt[nm][all][cn] );
       sten := Gen_Sten( evol_eq_LHS_tbl[nm][jj] -  tt[nm][all][cn]);
       if (sten <> 0 ) then
         Gen_Res_Code(sten,input="d",proc_name=LowerCase(cat("ire_",nm,cn)));
         Gen_Eval_Code(sten,input="d",proc_name=LowerCase(cat("irev_",nm,cn)));
         Other_Headers := Other_Headers , cat("ire_",nm,cn);
       end if;
    end do
  else
       all_expr:=all_expr + ( evol_eq_LHS_tbl[nm] -  tt[nm][all] );
       sten := Gen_Sten( evol_eq_LHS_tbl[nm] -  tt[nm][all]);
       if (sten <> 0) then
         Gen_Res_Code(sten,input="d",proc_name=LowerCase(cat("ire_",nm)));
         Gen_Eval_Code(sten,input="d",proc_name=LowerCase(cat("irev_",nm)));
         Other_Headers := Other_Headers , cat("ire_",nm);
       end if;
  end if;
end do:


text:
# Forward Kinematics for Robot based on MDH frames
# Introduction
# Direkte Kinematik basierend auf MDH-Parametern berechnen
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# floatb -> floating base wird durch base twist (Geschwindigkeit der Basis) oder vollständige Orientierung (Euler-Winkel) berücksichtigt
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# mdh_kinematics -> Berechnung der Vorwärtskinematik mit modifizierten DH-Parametern nach [KhalilKle1986]
#
# Prinzip:
# Berechne die direkte Kinematik.
# Zusätzlich können kinematische Zwangsbedingungen direkt berücksichtigt werden. Diese sorgen dafür, dass MDH-Winkel durch analytische Ausdrücke verallgemeinerter Koordinaten ersetzt werden.
# Authors
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-03
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Sources
# [KhalilKle1986] Khalil, W. & Kleinfinger, J.: A new geometric notation for open and closed-loop robots (1986)
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
codegen_act := true:
codegen_debug := false:
codegen_opt := 2:
# Code-Optimierung für die Kinematik (zusammenfassung paralleler Achsen).
# Ist für einige Systeme vorteilhaft (serielle Ketten mit parallelen Achsen), für andere nicht (z.B. Baumstrukturen mit parallelen Achsen. Hier kann die Code-Optimierung mehr aus nicht trigonometrisch optimierten Termen herausholen).
# Die Terme werden aber auf jeden Fall übersichtlicher.
codegen_kinematics_opt := true:
# Substitutionsreihenfolge von kinematischen Zwangsbedingungen (falls vorhanden) einstellen.
# 1: Ersetzung bereits in Einzel-Transformationsmatrizen (dann ist die Zusammenfassung für mehrere Achsen aus codegen_kinematics_opt nicht mehr so leicht möglich). Muss 1 sein für trigonometrische Ersetzung.
# 2: Ersetzung erst nach der Zusammenfassung der Einzel-Transformationsmatrizen zu (kumulierten) Gesamt-Transformationsmatrizen
codegen_kinematics_subsorder := 1:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../helper/proc_combine2":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
# Funktionen aus Robotik-Repo
read("../robotics_repo_path"):
read(sprintf("%s/transformation/maple/proc_eulxyz2r", robotics_repo_path)):
# Homogene Transformationsmatrix für Rotation aus XYZ-Euler-Winkel
eulxyztr := proc (alpha, beta, gamma)
  local T;
  T := Matrix(4,4,`<,>`(`<|>`(eulxyz2r(<alpha; beta; gamma>),`<,>`(0,0,0)),`<|>`(0,0,0,1)));
  return T
end proc:
read "../robot_codegen_definitions/robot_env":
printf("%s. Generiere Kinematik für %s\n", FormatTime("%Y-%m-%d %H:%M:%S"), robot_name):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
printf("%s. Alle Daten geladen. Generiere Kinematik für %s\n", FormatTime("%Y-%m-%d %H:%M:%S"), robot_name):
# Term-Vereinfachungen einstellen
if not assigned(simplify_options) or simplify_options(2)=-1 then
  use_simplify := 1: # standardmäßig simplify-Befehle anwenden. Sind meistens sinnvoll und in Ausnahmen nicht allzu schädlich
else
  use_simplify := simplify_options(2): # zweiter Eintrag ist für Kinematik
end if:

# Kinematische Zwangsbedingungen
# Lade Ausdrücke für kinematische Zwangsbedingungen (Verknüpfung von MDH-Gelenkwinkeln durch verallgemeinerte Koordinaten)
# Lese Variablen: kintmp_subsexp, kintmp_qs, kintmp_qt
kin_constraints_exist := false:
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
constrfile := sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
if FileTools[Exists](constrfile) then
  read constrfile:
end if:
if kin_constraints_exist = true then:
  kintmp_qs := kintmp_qs: # gelesene Variable sonst nicht sichtbar
  kintmp_qt := kintmp_qt: # gelesene Variable sonst nicht sichtbar
  kintmp_subsexp := kintmp_subsexp:
  printf("%s. Kinematische Zwangsbedingungen gelesen.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
else
  kintmp_qs := Matrix(1,1):
  kintmp_qt := Matrix(1,1):
  kintmp_subsexp := Matrix(1,2):
  kin_constraints_exist := false:
  # Zwangsbedingungen neu speichern, damit diese auch für andere Skripte verfügbar sind (als dummy-Variablen).
  save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
end if:
# Calculate Forward Kinematics (Single-Joint Transformation)
# Trf is the Matrix of Transformation from i-1 to i
# Trf_c is the cummulated Matrix of Transformation from 0 to i
printf("%s. Beginne Berechnung der Kinematik für einzelne Gelenk-Transformationen.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
Trf := Matrix(4, 4, NJ): # Diese Initialisierung bringt nichts (initialisiert nur 4x4-Matrix)
;
for i from 1 to NJ do
  Trf(1 .. 4, 1 .. 4, i) := Matrix(4,4):
end do:
for i from 1 to NJ do
  Trf(1 .. 4, 1 .. 4, i) := trafo_mdh_full(alpha[i,1], a[i,1], theta[i,1], d[i,1], beta[i,1], b[i,1]):
end do:
# Kinematische Zwangsbedingungen ersetzen (vor Berechnung der Gesamt-Transformation)
# Substituiere allgemeine Ausdrücke der Winkel der Parallelstruktur mit kinematischen Zwangsbedingungen in Abhängigkeit der Haupt-Gelenkwinkel
if kin_constraints_exist and codegen_kinematics_subsorder = 1 then
  for i from 1 to NJ do # Index über Transformationsmatrizen aller Körper
    for ix from 1 to 4 do # Index über Zeilen der Transformationsmatrizen
      for iy from 1 to 4 do # Index über Spalten der Transformationsmatrizen
        # Substituiere Zeitvariablen
        Trf(ix, iy, i) := convert_t_s( Trf(ix, iy, i) ):
        # Substituiere sin und cos der Winkel (einfachere Ausdrücke)
        Trf(ix, iy, i) := subs_kintmp_exp(Trf(ix, iy, i)):
        # Substituiere die verbleibenden Winkel direkt (einige Winkel sind nicht in den Ersetzungsausdrücken enthalten, da sie keine problematischen arctan-Ausdrücke enthalten.
        for jj from 1 to RowDimension(kintmp_qt) do # Index über zu ersetzende Winkel
          Trf(ix, iy, i) := subs( { kintmp_s(jj, 1) = kintmp_qs(jj, 1) }, Trf(ix, iy, i) ):
        end do:
        Trf(ix, iy, i) := convert_s_t( Trf(ix, iy, i) ):
      end do:
    end do:
    if use_simplify>=1 then
    	 tmp_l1 := length(Trf(1 .. 4, 1 .. 4, i)): tmp_t1:=time():
    	 Trf(1 .. 4, 1 .. 4, i) := simplify2(Trf(1 .. 4, 1 .. 4, i)):
    	 tmp_l2 := length(Trf(1 .. 4, 1 .. 4, i)): tmp_t2:=time():
    	 printf("%s: Term für (Einzel-)Trafo-Matrix %d->%d vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    	   FormatTime("%Y-%m-%d %H:%M:%S"), v(i), i, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
    end if:
  end do:
  printf("%s. Ersetzungen der MDH-Parameter mit Ergebnissen der Parallelstruktur in verallgemeinerten Koordinaten erfolgreich (vor Berechnung der Gesamt-Transformation).\n", FormatTime("%Y-%m-%d %H:%M:%S")):
end if:
# Calculate Forward Kinematics (Multi-Joint Transformation)
Trf_c := Matrix(4, 4, NJ+1): # Diese Initialisierung bringt nichts (initialisiert nur 4x4-Matrix)
for i from 1 to NJ+1 do
  Trf_c(1 .. 4, 1 .. 4, i) := Matrix(4,4): # Vollständige Initialisierung
end do:
# Basis-Transformation: Unterschiedliche Darstellungsmethoden. Führen zu unterschiedlichen verallgemeinerten Koordinaten.
if base_method_name = "twist" then:
  Trf_c(1 .. 4, 1 .. 4, 1) := transl(X_base_t[1..3,1]):
end:
if base_method_name = "eulxyz" then:
  Trf_c(1 .. 4, 1 .. 4, 1) := transl(X_base_t[1..3,1]) . eulxyztr(X_base_t[4,1], X_base_t[5,1], X_base_t[6,1]):
end:
printf("%s. Beginne direkte Kinematik. Nutze die Methode %s für die Basis-Orientierung.\n", FormatTime("%Y-%m-%d %H:%M:%S"), base_method_name):
# Kinematik aller Körper mit MDH-Ansatz Bestimmen. [KhalilKle1986].
# Einfache Kinematik: Multiplikation an vorherige Transformationsmatrix
if not(codegen_kinematics_opt) then
  printf("%s. Berechne Vorwärts-Transformation der Kinematik ohne spezielle Optimierung.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  for i from 1 to NJ do
    # Index des vorherigen Koordinatensystems
    j := v(i)+1:
    Trf_c(1 .. 4, 1 .. 4, i+1) := Multiply(Trf_c(1 .. 4, 1 .. 4, j), Trf(1 .. 4, 1 .. 4, i)):
    if use_simplify>=1 then
    	 tmp_l1 := length(Trf_c(1 .. 4, 1 .. 4, i+1)): tmp_t1:=time():
    	 Trf_c(1 .. 4, 1 .. 4, i+1) := simplify2(Trf_c(1 .. 4, 1 .. 4, i+1)):
    	 tmp_l2 := length(Trf_c(1 .. 4, 1 .. 4, i+1)): tmp_t2:=time():
    	 printf("%s: Term für (Gesamt-)Trafo-Matrix 0->%d vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    	   FormatTime("%Y-%m-%d %H:%M:%S"), i, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
    end if:
  end do:
end if:
#Trf_c_orig := copy(Trf_c): # Debug-Ausdruck zum Vergleich der obigen mit der unteren Methode
;
# Optimierung der Terme: Summentheorem für Drehungen um aufeinanderfolgende Achsen.
# Bei Verwendung dieser Optimierung funktioniert die Substitution mit trigonometrischen Ausdrücken von abhängigen Gelenken eventuell anders.
# Die Funktionsweise kann mit aktivierten print-Befehlen nachvollzogen werden.

if codegen_kinematics_opt then
  printf("%s. Beginne mit der Optimierung der Kinematik für parallele Gelenke.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  for i from 1 to NJ do
    #printf("----------------------\ni=%d\n",i):
    # Bestimme Vorwärtskinematik für das Koordinatensystem i (0=Basis)
    # Gehe solande rückwärts, wie die Gelenkachsen parallel sind
    j := i: # Lauf-Index über Vorgänger-Elemente (0=Basis)
    j2 := 0:# Index des Segmentes (0=Basis), das vor der Kette von parallelen Achsen ist
    Trf_tmp := IdentityMatrix(4,4); # Transformation von diesem Segment zum aktuellen
    Kette_akt := []:
    #printf("Trf_tmp neu initialisiert:\n"):
    #print(Trf_tmp);
    for k from 1 to NJ do # Schleife mit Dummy-Länge. Wird abgebrochen, falls beendet.
      #printf("j=%d\n",j):
      Kette_akt := [j, op(Kette_akt)]:
      # Additionstheorem für Drehung um parallele Achsen (combine)
      # Prüfe, welche Form rechentechnisch am günstigsten ist und wähle diese (wrapper combine2)
      Trf_tmp := simplify2(combine2(Matrix(Trf(1 .. 4, 1 .. 4, j) . Trf_tmp))):
      #printf("Trf_tmp aktualisiert mit Trf %d. Aktuelle Kette: %s\n", j, convert(Kette_akt, string)):
      #print(Trf_tmp);
      if not(alpha(j) = 0) then
        # Die vorherige Achse ist nicht parallel zu dieser
        # Weitere Vereinfachungen ergeben keinen Sinn
        j2 := v(j):
        #printf("alpha ungleich null. Abbruch j2=%d.\n", j2):
        break:
      #else
      #  printf("alpha ist Null. Vorherige Achse ist parallel\n"):
      end if:
      if v(j) = 0 then
        break: # Vorgänger ist Basis. Keine weitere Untersuchung
      end if:
      #printf("v(%d)=%d\n", j, v(j)):
      j := v(j): # Nehme das Vorgänger-Segment
    end do:
    Trf_c(1 .. 4, 1 .. 4, i+1) := Matrix(Trf_c(1 .. 4, 1 .. 4, j2+1)) . Trf_tmp:
    if use_simplify>=1 then
    	 tmp_l1 := length(Trf_c(1 .. 4, 1 .. 4, i+1)): tmp_t1:=time():
    	 Trf_c(1 .. 4, 1 .. 4, i+1) := simplify2(Trf_c(1 .. 4, 1 .. 4, i+1)):
    	 tmp_l2 := length(Trf_c(1 .. 4, 1 .. 4, i+1)): tmp_t2:=time():
    	 printf("%s: Term für (Gesamt-)Trafo-Matrix 0->%d vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs\n", \
    	   FormatTime("%Y-%m-%d %H:%M:%S"), i, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
    end if:
    #printf("Trf_tmp an Trf_c angehängt (Eintrag zu Körper %d. %d -> %d mit Kette %s)\n", i, j2, i, convert(Kette_akt, string)):
    #print(Trf_c(1 .. 4, 1 .. 4, i+1)):
  end do:
  printf("%s. Berechnung der direkten Kinematik mit Vereinfachungen für parallele Achsen erfolgreich.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S")
):
end if:
# Kinematische Zwangsbedingungen ersetzen (nach Berechnung der Gesamt-Transformation)
# Substituiere allgemeine Ausdrücke der Winkel der Parallelstruktur mit kinematischen Zwangsbedingungen in Abhängigkeit der Haupt-Gelenkwinkel
# Da die Gesamt-Transformation Trf_c an dieser Stelle schon berechnet wurde, müssen die Ausdrücke hier in Trf und Trf_c ersetzt werden.
if kin_constraints_exist and codegen_kinematics_subsorder = 2 then
  printf("%s. Beginne Einsetzung der kinematischen Zwangsbedingungen in die direkte Kinematik (Gesamt-Transformation).\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  for i from 1 to NJ do # Index über Transformationsmatrizen aller Körper
    for ix from 1 to 4 do # Index über Zeilen der Transformationsmatrizen
      for iy from 1 to 4 do # Index über Spalten der Transformationsmatrizen
        # Substituiere Zeitvariablen
        Trf(  ix, iy, i)   := convert_t_s( Trf(  ix, iy, i) ):
        Trf_c(ix, iy, i+1) := convert_t_s( Trf_c(ix, iy, i+1) ):
        # Substituiere sin und cos der Winkel (einfachere Ausdrücke)
        Trf(ix,   iy, i)   := subs_kintmp_exp(Trf(  ix, iy, i)):
        Trf_c(ix, iy, i+1) := subs_kintmp_exp(Trf_c(ix, iy, i+1)):
        # Substituiere die verbleibenden Winkel direkt (einige Winkel sind nicht in den Ersetzungsausdrücken enthalten, da sie keine problematischen arctan-Ausdrücke enthalten.
        for jj from 1 to RowDimension(kintmp_qt) do # Index über zu ersetzende Winkel
          Trf(ix,   iy, i)   := subs( { kintmp_s(jj, 1) = kintmp_qs(jj, 1) }, Trf(  ix, iy, i) ):
          Trf_c(ix, iy, i+1) := subs( { kintmp_s(jj, 1) = kintmp_qs(jj, 1) }, Trf_c(ix, iy, i+1) ):
        end do:
        Trf(ix,   iy, i)   := convert_s_t( Trf(  ix, iy, i) ):
        Trf_c(ix, iy, i+1) := convert_s_t( Trf_c(ix, iy, i+1) ):
      end do:
    end do:
    if use_simplify>=1 then
    	 tmp_l1 := length(Trf(1 .. 4, 1 .. 4, i)): tmp_t1:=time():
    	 Trf(1 .. 4, 1 .. 4, i) := simplify2(Trf(1 .. 4, 1 .. 4, i)):
    	 tmp_l2 := length(Trf(1 .. 4, 1 .. 4, i)): tmp_t2:=time():
    	 printf("%s. Term für Einzel-Trafo-Matrix %d->%d vereinfacht. Länge: %d->%d Rechenzeit %1.1fs.\n", \
    	   FormatTime("%Y-%m-%d %H:%M:%S"), v(i), i, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
    	 tmp_l1 := length(Trf_c(1 .. 4, 1 .. 4, i+1)): tmp_t1:=time():
    	 Trf_c(1 .. 4, 1 .. 4, i+1) := simplify2(Trf_c(1 .. 4, 1 .. 4, i+1)):
    	 tmp_l2 := length(Trf_c(1 .. 4, 1 .. 4, i+1)): tmp_t2:=time():
    	 printf("%s. Term für Gesamt-Trafo-Matrix 0->%d vereinfacht. Länge: %d->%d Rechenzeit %1.1fs.\n", \
    	   FormatTime("%Y-%m-%d %H:%M:%S"), i, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
    end if:
  end do:
  printf("%s. Ersetzungen der MDH-Parameter mit Ergebnissen der Parallelstruktur in verallgemeinerten Koordinaten erfolgreich (nach Berechnung der Gesamt-Transformation).\n", \
    FormatTime("%Y-%m-%d %H:%M:%S")):
end if:
# Export
printf("%s. Berechnung der Kinematik beendet. Beginne Matlab-Export.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
# Maple-Export
save Trf, Trf_c, sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):

# Export des symbolischen Ausdrucks für alle kumulierten Transformationsmatrizen auf einmal.
# Export ohne letzte Zeile mit [0 0 0 1]
Trf_c_Export := Matrix((NJ+1)*3, 4):
for i from 1 to NJ+1 do
  Trf_c_Export((i-1)*3+1 .. 3*i, 1..4) := Trf_c(1..3, 1..4, i):
end do:
if codegen_act then
  MatlabExport(convert_t_s(Trf_c_Export), sprintf("../codeexport/%s/tmp/fkine_mdh_floatb_%s_rotmat_matlab.m", robot_name, base_method_name), codegen_opt):
end if:

# Export des symbolischen Ausdrucks für alle Gelenk-Transformationsmatrizen auf einmal.
# Export ohne letzte Zeile mit [0 0 0 1]
Trf_Export := Matrix((NJ)*3, 4):
for i from 1 to NJ do
  Trf_Export((i-1)*3+1 .. 3*i, 1..4) := Trf(1..3, 1..4, i):
end do:
if codegen_act then
  MatlabExport(convert_t_s(Trf_Export), sprintf("../codeexport/%s/tmp/joint_transformation_mdh_rotmat_matlab.m", robot_name), codegen_opt):
end if:

# Export des symbolischen Ausdrucks für jede Transformationsmatrix einzeln
for i from 1 to NJ+1 do
  if codegen_act then
    MatlabExport(convert_t_s(Trf_c(1 .. 4, 1 .. 4, i)), sprintf("../codeexport/%s/tmp/fkine_%d_floatb_%s_rotmat_matlab.m", robot_name, i-1, base_method_name), codegen_opt):
  end if:
end do:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


text:
#@ Not autoload

# This file relates to the cubic with equation x[1]^3 + x[2]^3 + x[3]^3 = nu * x[1] * x[2] * x[3]

om := exp( 2*Pi*I/3);
ob := exp(-2*Pi*I/3);
R := [1,om,ob];

# The group G acts on affine 3-space, preserving the above cubic.
# The quotient G/Z acts on the projective plane, preserving our cubic curve

id := IdentityMatrix(3);
rho := <<0|0|1>,<1|0|0>,<0|1|0>>:
zeta := DiagonalMatrix([1,om,ob]):

Z := [id,om * id, ob * id];
T := [seq(seq(DiagonalMatrix(rationalize(expand([om^i,om^j,om^(-i-j)]))),j=0..2),i=0..2)]:
S := map(perm_matrix,combinat[permute](3)):
G := [seq(seq(t.s,t in T),s in S)]:

# rho and zeta commute in G/Z
for i from -1 to 1 do
 for j from -1 to 1 do
  g[i,j] := map(rationalize,rho^i . zeta^j);
 od;
od:

# We now define matrices U[i,j] = <e[i,j] | u[i,j] | v[i,j]>
# Each e[i,j] generates a 1-dimensional summand in affine 3-space
# On the complement of the locus where nu = 3 = 0, each pair
U[0,0] := <<1|1|nu>,<-1|-nu-1|3>,<0|3|nu>>;
for i from -1 to 1 do
 for j from -1 to 1 do
  U[i,j] := map(expand,g[i,j] . U[0,0]);
  e[i,j] := convert(Column(U[i,j],1),list);
  u[i,j] := convert(Column(U[i,j],2),list);
  v[i,j] := convert(Column(U[i,j],3),list);
 od:
od:


text:
f1(){
	[x, y] = stack(), stack();
	x*7-y; // return value
}
f2()inline{
	[x, y] = stack(), stack();
	x*7-y; // return value
}
f3()inline, CRPL{swap 7 mul swap sub}
f4()CRPL{swap 7 mul swap sub}

main()inline{
	ShowTraceLog();
	ClearTraceLog();

	Trace("f1 " $ f1(1, 2));
	Trace("f2 " $ f2(1, 2));
	Trace("f3 " $ f3(1, 2));
	Trace("f4 " $ f4(1, 2));
}