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_c *)

$include "gga_c_gapc.mpl"

(* override definition of gap_C *)
gap_G := (rs, z, xt, par) -> RS_FACTOR^2/8 * xt^2/rs^2:

f_kcis := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  + f_gap(rs, z, xt)
  - xs0^2/(8*ts0) * (1 + z)/2 * f_gap(rs,  1, xs0)
  - xs1^2/(8*ts1) * (1 - z)/2 * f_gap(rs, -1, xs1):

f  := (rs, z, xt, xs0, xs1, ts0, ts1, us0, us1) ->
  f_kcis(rs, z, xt, xs0, xs1, ts0, ts1):


text:
program assignchar;
var x, y : char;

ext: mpl

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:
(*
 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: mgga_exc *)

cs_a := -0.04918:
cs_b :=  0.132:
cs_c :=  0.2533/RS_FACTOR:
cs_d :=  0.349/RS_FACTOR:

cs_thf := (z, u, t) ->
  opz_pow_n(z,8/3)*2^(-8/3)*(t - u/8):

(* This is Equation (15) of Lee1988_785 *)
(* Note that gamma = (1 - z^2) *)
f_cs := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) ->
  cs_a*(1 - z^2)/(1 + cs_d*rs) * (1 + 2*cs_b*exp(-cs_c*rs)*(
    cs_thf(z, u0, t0) + cs_thf(-z, u1, t1) - t_vw(z, xt, u0, u1)
  )):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) ->
  f_cs(rs, z, xt, xs0, xs1, u0, u1, t0, t1):


text:
# This file is about Singer systems of endomorphisms of the additive formal
# group.  Such a system has groups G(i,j) (for 0 <= i <= j <= n say),
# all of which are equal to the standard additive formal group.  There
# are morphisms G(j,k) -> G(i,l) whenever i <= j <= k <= l, with the
# evident functoriality condition.  The horizontal maps are linear
# and the vertical maps are short isogenies.  Every subdiagram
# G(i,i) -> G(i,j) <- G(j,j) is of the form G -> G/a <- G//a.

W0V := (k) -> mul(V[i],i=0..k-1):

WW := proc(k,l,m)
  local J,P;
  J := [seq(j,j=k..l-1)];
  P := combinat[choose](J,l-k-m);
  return add(mul(V[u[t]]^(2^(k+m+t-1-u[t])),t=1..l-k-m),u in P);
end:

WWW := (i,j,k,l,m) -> WW(k,l,m) * W0V(j)^(2^(k+m-j)) / W0V(i)^(2^(l-i)):

# G(j,k) to G(i,l)
f := proc(i,j,k,l)
 local x;
 return unapply(sort(collect(expand(add(WWW(i,j,k,l,m) * x^(2^m),m=0..l-k)),x)),x):
end:

vV_rule := {seq(v[i]=V[i]/W0V(i),i=0..20)}:

Vv_rule := {seq(
 V[i] = v[i] * mul(v[j]^(2^(i-j-1)),j=0..i-1),
 i=0..20)};

check_id := proc(i::nonnegint,k::nonnegint)
 local x;

 if not (i <= k) then return FAIL; fi;

 _ASSERT(f(i,i,k,k)(x) = x,"identity maps",[i,k]);

end:

check_compose := proc(h::nonnegint,i::nonnegint,j::nonnegint,
                      k::nonnegint,l::nonnegint,m::nonnegint)
 local x,err;

 if not(h <= i and i <= j and j <= k and k <= l and l <= m) then
  return FAIL;
 fi;

 err := modp(expand(f(h,i,l,m)(f(i,j,k,l)(x)) - f(h,j,k,m)(x)),2);

 _ASSERT(err = 0,"maps compose correctly",[h,i,j,k,l]);
end:

check_linear := proc(i::nonnegint,j::nonnegint,k::nonnegint)
 local x,y,err;

 if not(i <= j and j <= k) then
  return FAIL;
 fi;

 y := f(i,j,k,k)(x);
 err := expand(y - coeff(y,x,1) * x);

 _ASSERT(err = 0,"horizontal map is linear",[i,j,k]);
end:

check_short_isogeny := proc(i::nonnegint,k::nonnegint,l::nonnegint)
 local x,y,err;

 if not(i <= k and k <= l) then
  return FAIL;
 fi;

 y := f(i,i,k,l)(x);

 _ASSERT(degree(y,x) = 2^(l-k) and coeff(y,x,2^(l-k)) = 1,
         "vertical map is a short isogeny",[i,k,l]);
end:

check_dilation := proc(i::nonnegint,p::nonnegint)
 local x,y,z,err;

 y := f(i,i,i,i+p)(x);
 z := f(i,i+p,i+p,i+p)(x);
 _ASSERT(expand(coeff(y,x,1) - coeff(z,x,1)) = 0,"dilation condition",[i,j,k]);
end:

check_vV_rules := proc(n)
 _ASSERT(
  `and`(seq(evalb(subs(vV_rule,subs(Vv_rule,V[i])) = V[i]),i=0..n)) and
  `and`(seq(evalb(subs(Vv_rule,subs(vV_rule,v[i])) = v[i]),i=0..n)),
  "rules for converting v[i] to V[i] and vice versa"
 );

end:


text:
read "print_comp.mpl":
all_var_set:=[
   [DPH,0,dn] , [dg,2,dn] ,
         [met,2,dn], [invmet,2,up], [DDL,2,dn],
   [DTK,0,dn],
         [RR,2,dn], [U,2,dn] , [P,2,dn],
   [DAA,2,dn],
   [DLAM,1,up],
   [C1,0,dn],
   [C5,0,dn],
   [HS,0,dn]
             ];
s:="";
for i from 1 to nops(all_var_set) do
  s:=cat(s,print_comp(all_var_set[i][1],[x,theta,phi],all_var_set[i][2],diag=true,tp=all_var_set[i][3])):
end do:
FP:=fopen("allequations.mpl",WRITE);
writeline(FP,s);
fclose(FP);


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 *)

coeffs_beefvdw := [
   1.516501714e0,   4.413532099e-1, -9.182135241e-2, -2.352754331e-2,
   3.418828455e-2,  2.411870076e-3, -1.416381352e-2,  6.975895581e-4,
   9.859205137e-3, -6.737855051e-3, -1.573330824e-3,  5.036146253e-3,
  -2.569472453e-3, -9.874953976e-4,  2.033722895e-3, -8.018718848e-4,
  -6.688078723e-4,  1.030936331e-3, -3.673838660e-4, -4.213635394e-4,
   5.761607992e-4, -8.346503735e-5, -4.458447585e-4,  4.601290092e-4,
  -5.231775398e-6, -4.239570471e-4,  3.750190679e-4,  2.114938125e-5,
  -1.904911565e-4,  7.384362421e-5
]:

k  := 4:
xi := p -> 2*p/(k + p) - 1:

with(orthopoly):
f  := x -> add(coeffs_beefvdw[i]*P(i-1, xi(X2S^2*x^2)), i=1..30):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
(*
 Copyright (C) 2020 Susi Lehtola

 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: mgga_exc *)

gea2_s := x -> X2S*x:
gea2_q := u -> X2S^2*u:

gea2_f0 := (s, q) -> 1 + 5/27*s^2 + 20/9*q:
gea2_f := (x, u) ->
  gea2_f0(gea2_s(x), gea2_q(u)):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) ->
  mgga_kinetic(gea2_f, rs, z, xs0, xs1, u0, u1):


text:
# Berechne kinematische Zwangsbedingungen für den Palettierroboter KUKA KR700PA aus 2 Viergelenkketten
# Einleitung
# Die kinematischen Zwangsbedingungen werden als Ersatzungsausdruck für die abhängigen Winkel aufgestellt.
#
# palh3m1TE->KUKA KR700PA, modellierung der Zwangsbedingungen mit Ausdrücken für trigonometrische Elimination
# kinematic_constraint->Kinematische Zwangsbedingungen
# Datenblatt des Roboters:
# unter :   https://www.kuka.com/en-de/products/robot-systems/industrial-robots/kr-700-pa
# Voraussetzung(sehr wichtig!!!!!):
# Viergelenkkette muss mit der Toolbox berechnet worden sein (Arbeitsblatt "fourbar1TE_kinematic_constrains.mv")
# Quelle:
# SA Shan: Shan2019_S828, "Reduktion der Modellkomplexität von seriell-hybriden Palettierrobotern"
# Autor:
# Weipu Shan
# Studienarbeit bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2018-08
# (C) Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
# Import (Library)
interface(warnlevel=0): 		# Unterdrücke die folgende Warnung.
restart: 				# Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true:               # für Speicherung
;
with(StringTools):                                  # für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
#with(ListTools):
codegen_act := true:
codegen_opt := 1:      # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0:    # Zur Code-Generierung auch für Nicht-Inert Ausdrücke
;
# Import(hybriddyn)
read "../helper/proc_MatlabExport":    #Matlab Export
read "../transformation/proc_rotx":      #Rotation um X
read "../transformation/proc_roty":      #Rotation um Y
read "../transformation/proc_rotz":      #Rotation um Z
read "../helper/proc_convert_s_t":        #zeitabhängige Variable für Ableitung
read "../helper/proc_convert_t_s":        #konstante Variable
read "../robot_codegen_constraints/proc_subs_kintmp_exp": #substitute Spalte 1 mit Spalte 2
with (RealDomain):            # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":                                                      #aktuelle Roboter, MDH-Tabelle
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions",robot_name):
# Variable Initialisierung
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten  (Nur Aktive Gelenke)
kintmp_qs := Matrix(RowDimension(kintmp_s),1):       	#kintmp_s von aktuellen Roboter, MDH
;
#
# Konstante Winkel bereits hineinschreiben (Passive Gelenk)
for i from 1 to RowDimension(kintmp_s) do
	if diff(kintmp_s(i,1),t) = 0 then
		#konstante Variable, passive Winkel von aktuellem Roboter, MDH
		kimtmp_qs(i,1) := kintmp_s(i,1):
  	end if:
end do
;
#
# Variable definieren für die Hilfswinkel
# Ersetzungsausdrücke definieren
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zusuchender Ausdruck (sin oder cos eines Winkels)	#第一列是要找到的三角函数表达
# Spalte 2: Einzusetzender Ausdruck			#第二列是代替的表达
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
# Ähnliches Vorgehen wie in [1].
#
# #新建kintmp_subsexp 12X2, 第一列被动轴是sin cos 形式，第二列初始化和第一列相同，当前机器人的sin cos
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do 			#kintmp_s von aktuellen Roboter
	kintmp_subsexp(2*i-1,1) := sin(kintmp_s(i,1)):
	kintmp_subsexp(2*i,  1) := cos(kintmp_s(i,1)):
	# Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
	kintmp_subsexp(2*i-1,2) := kintmp_subsexp(2*i-1, 1):
	kintmp_subsexp(2*i,  2) := kintmp_subsexp(2*i,   1):
end do:
#
# Sicherung der aktuelle Daten
backup_kintmp_qs := kintmp_qs:		#Backup von aktuellem Roboter, konstante Variable
backup_kintmp_qt := kintmp_qt:		#Backup von aktuellem Roboter, zeitabhängige Variable
backup_kintmp_subsexp := kintmp_subsexp:   	#Backup von aktuellem Roboter, sin, cos von 4 Winkeln

;
# Daten von Fourbar übernehmen
read sprintf("../codeexport/fourbar1TE/tmp/kinematic_constraints_maple_inert.m"):
# Winkeln von Fourbar会先覆盖原先的变量 kintmp_subsexp, kintmp_qs, kintmp_qt,所以一开始需要备份

;
# Speichern WInkeln von Fourbar, mit Fourbar Variable
kintmp_subsexp_fourbar := kintmp_subsexp:  	#sin,cos von 4 Winkeln, Fourbar
kintmp_qs_fourbar := kintmp_qs:		#konstante Variable, Fourbar
kintmp_qt_fourbar := kintmp_qt:		#zeitabhängige Variable, Fourbar
;
#
# 以下内容有疑问 Was ist der Sinn darunten
#kintmp_subsexp_fourbar:
#winkel(6,2):
#kintmp_subsexp_fourbar(6,2):
#kintmp_subsexp_fourbar(3,3):
#winkel_neu := Matrix(6,2):
winkel_neu(1,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(3,2)):
winkel_neu(2,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(4,2)):
winkel_neu(3,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(1,2)):
winkel_neu(4,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(2,2)):
winkel_neu(5,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(5,2)):
winkel_neu(6,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(6,2)):
#
# Vorherige Werte für dieses System wieder herstellen
kintmp_qs := backup_kintmp_qs:		#aktuellem Roboter, konstante Variable
kintmp_qt := backup_kintmp_qt:		#aktuellem Roboter, zeitabhängige Variable
kintmp_subsexp := backup_kintmp_subsexp:		#aktuellem Roboter, sin, cos von 4 Winkeln

;
winkel_neu(1..6,1);
# Viergelenkkette 1: A-B-C-D
# Von ABCD: eta1, eta3, eta4  -----> Von Fourbar: rho, eta, xi
winkel1 := <<sin_eta1|sin(eta1_s)>;<cos_eta1|cos(eta1_s)>;<sin_eta3|sin(eta3_s)>;<cos_eta3|cos(eta3_s)>;<sin_eta4|sin(eta4_s)>;<cos_eta4|cos(eta4_s)>>:
# Ausdrücke von Fourbar in aktuellem Roboter mitbringen
for i from 1 to 6 do
	winkel1(i,2) := kintmp_subsexp_fourbar(i,2):       		#Subsexpress von Fourbar übernehmen
	cos_phi_s := sin(qJ2s)*sin(phi2) - cos(qJ2s)*cos(phi2):	#phi_s ist unabhängige Winkel von VGK, hier ist phi_s = eta2_s
	sin_phi_s := sin(qJ2s)*cos(phi2) + cos(qJ2s)*sin(phi2):	#qJ2s ist echte aktive Winkel von Roboter, eta2 ^= phi = Pi - qJ2 - phi2
	winkel1(i,2) := subs({sin(qJ1s)=sin_phi_s},winkel1(i,2)):
	winkel1(i,2) := subs({cos(qJ1s)=cos_phi_s},winkel1(i,2)):
	winkel1(i,2) := subs({l1=AB}, winkel1(i,2)):
	winkel1(i,2) := subs({l2=DA}, winkel1(i,2)):
	winkel1(i,2) := subs({l3=DC}, winkel1(i,2)):
	winkel1(i,2) := subs({l4=BC}, winkel1(i,2)):
end do:
winkel1(1..6,1):
# eta1, eta3, eta4 wird wie rho, eta, xi in Fourbar darstellen
sin_eta1_s := winkel1(1,2):	#eta1 = qJ2s = rho
;
cos_eta1_s := winkel1(2,2):
sin_eta3_s := winkel1(3,2):	#eta3 = qJ3s = eta
;
cos_eta3_s := winkel1(4,2):
sin_eta4_s := winkel1(5,2):	#eta4 = qJ4s = xi
;
cos_eta4_s := winkel1(6,2):
# eta16, eta79, eta27, abhängige Winkeln von aktuellem Roboter
cos_eta16_s := cos_eta1_s*cos(phi1) + sin_eta1_s*sin(phi1):	#eta16 = eta1 - phi1
;
sin_eta16_s := sin_eta1_s*cos(phi1) - cos_eta1_s*sin(phi1):
cos_eta79_s := - cos_eta4_s:				#eta79 = Pi-eta4
;
sin_eta79_s := sin_eta4_s:
cos_eta27_s := cos_eta3_s*cos(phi79) + sin_eta3_s*sin(phi79):	#eta27 = eta3 - phi79
;
sin_eta27_s := sin_eta3_s*cos(phi79) - cos_eta3_s*sin(phi79):
# Viergelenkkette 2: B-E-F-G
winkel2 := <<sin_xi1|sin(xi1_s)>;<cos_xi1|cos(xi1_s)>;<sin_xi3|sin(xi3_s)>;<cos_xi3|cos(xi3_s)>;
<sin_xi4|sin(xi4_s)>;<cos_xi4|cos(xi4_s)>>:
# Ausdrücke von Fourbar in aktuellem Roboter mitbringen
for i from 1 to 6 do
	winkel2(i,2) := kintmp_subsexp_fourbar(i,2):       			#Subsexpress von Fourbar übernehmen
	#phi_s ist unabhängige Winkel von VGK, hier ist phi_s = xi2_s
	#qJ3s ist echte aktive Winkel von Roboter, xi2 ^= phi = Pi - qJ3 - (phi78-eta27)
	cos_phi_s := -cos(phi78+phi79)*(cos_eta3_s*cos(qJ3s)+sin_eta3_s*sin(qJ3s))-sin(phi78+phi79)*(sin_eta3_s*cos(qJ3s)-cos_eta3_s*sin(qJ3s)):
	sin_phi_s := sin(phi78+phi79)*(cos_eta3_s*cos(qJ3s)+sin_eta3_s*sin(qJ3s))-cos(phi78+phi79)*(sin_eta3_s*cos(qJ3s)-cos_eta3_s*sin(qJ3s)):
	#
	winkel2(i,2) := subs({sin(qJ1s)=sin_phi_s},winkel2(i,2)):
	winkel2(i,2) := subs({cos(qJ1s)=cos_phi_s},winkel2(i,2)):
	winkel2(i,2) := subs({l1=BG}, winkel2(i,2)):
	winkel2(i,2) := subs({l2=BE}, winkel2(i,2)):
	winkel2(i,2) := subs({l3=EP}, winkel2(i,2)):
	winkel2(i,2) := subs({l4=GP}, winkel2(i,2)):
end do:
winkel2(1..6,1):
# xi1, xi3, xi4 wird wie rho, eta, xi in Fourbar darstellen
sin_xi1_s := winkel2(1,2):	#xi1 = qJ2s = rho
;
cos_xi1_s := winkel2(2,2):
sin_xi3_s := winkel2(3,2):	#xi3 = qJ3s = eta
;
cos_xi3_s := winkel2(4,2):
sin_xi4_s := winkel2(5,2):	#xi4 = qJ4s = xi
;
cos_xi4_s := winkel2(6,2):
# xi78, xi34, xi410, abhängige Winkeln von aktuellem Roboter
cos_xi78_s := cos_xi1_s:				#xi78 = xi1
;
sin_xi78_s := sin_xi1_s:
cos_xi34_s := cos_xi3_s*cos(phi410) + sin_xi3_s*sin(phi410):	#xi34 = xi3 - phi410
;
sin_xi34_s := sin_xi3_s*cos(phi410) - cos_xi3_s*sin(phi410):
cos_xi410_s := -cos_xi4_s:				#xi410 = Pi - xi4
;
sin_xi410_s := sin_xi4_s:
# Kintmp_subsexp
kintmp_subsexp(1, 2) := sin_xi34_s:
kintmp_subsexp(2, 2) := cos_xi34_s:
kintmp_subsexp(3, 2) := sin_eta16_s:
kintmp_subsexp(4, 2) := cos_eta16_s:
kintmp_subsexp(5, 2) := sin_eta27_s:
kintmp_subsexp(6, 2) := cos_eta27_s:
kintmp_subsexp(7, 2) := sin_xi78_s:
kintmp_subsexp(8, 2) := cos_xi78_s:
kintmp_subsexp(9, 2) := sin_eta79_s:
kintmp_subsexp(10, 2) := cos_eta79_s:
kintmp_subsexp(11, 2) := sin_xi410_s:
kintmp_subsexp(12, 2) := cos_xi410_s:
kintmp_subsexp:

# Kintmp_qs
kintmp_qs(1,1) := %arctan(sin_xi34_s,cos_xi34_s):
kintmp_qs(2,1) := %arctan(sin_eta16_s,cos_eta16_s):
kintmp_qs(3,1) := %arctan(sin_eta27_s,cos_eta27_s):
kintmp_qs(4,1) := %arctan(sin_xi78_s,cos_xi78_s):
kintmp_qs(5,1) := %arctan(sin_eta79_s,cos_eta79_s):
kintmp_qs(6,1) := %arctan(sin_xi410_s,cos_xi410_s):

kintmp_subsexp:
# Export
winkel := <<sin_xi34|sin_xi34_s>; <cos_xi34|cos_xi34_s>; <sin_eta16 | sin_eta16_s>; <cos_eta16 | cos_eta16_s>;<sin_eta27 | sin_eta27_s>; <cos_eta27 | cos_eta27_s>;<sin_xi78|sin_xi78_s>;<cos_xi78|cos_xi78_s>;<sin_eta79|sin_eta79_s>;<cos_eta79|cos_eta79_s>;<sin_xi410|sin_xi410_s>;<cos_xi410|cos_xi410_s>>:
save winkel, sprintf("../codeexport/palh3m1TE/tmp/kinematic_constraints_maple_inert"):

kintmp_qt := convert_s_t(kintmp_qs):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name)
;
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name)
;
#printf("Ausdrücke für kintmp_subsexp gespeichert (Maple). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d",robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert" ,robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
#printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
interface(rtablesize=20):
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp ));
#kc_symbols :=Transpose(kc_symbols);
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
kc_symbols:
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);
#printf("Fertig. %s. CPU-Zeit bis hier: %1.2fs. \n", FormatTime("%Y-%m-%d %H:%M:%S"), time() - st):

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


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 *)

$define gga_c_pbe_params
$include "gga_c_pbe.mpl"

$include "lda_c_2d_amgb.mpl"
$undef xc_dimensions_2d

rs2D_factor := 1.704:
q2d_dd := 1e6:

q2d_rs2D := (rs, xt) -> rs2D_factor*rs*sqrt(X2S*xt)/RS_FACTOR:

q2d_fac := t -> t^4*(1 + t^2)/(q2d_dd + t^6):

q2d_f := (rs, z, xt, xs0, xs1) ->
  (1 - q2d_fac(tt(rs, z, xt)))*f_pbe(rs, z, xt, xs0, xs1) + q2d_fac(tt(rs, z, xt))*f_amgb(q2d_rs2D(rs, xt), z):

f := (rs, z, xt, xs0, xs1) -> q2d_f(rs, z, xt, 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: gga_exc *)
(* prefix:
  gga_x_pw91_params *params;

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

$ifdef gga_x_pw91_params
params_a_a     :=   0.19645:
params_a_b     :=   7.7956:
params_a_c     :=   0.2743:
params_a_d     :=  -0.1508:
params_a_f     :=   0.004:
params_a_alpha := 100:
params_a_expo  :=   4:
$endif

pw91_num := s -> (params_a_c + params_a_d*exp(-params_a_alpha*s^2))*s^2
         - params_a_f*s^params_a_expo:
pw91_den := s -> 1 + s*params_a_a*arcsinh(params_a_b*s) + params_a_f*s^params_a_expo:

pw91_f  := x -> 1 + pw91_num(X2S*x)/pw91_den(X2S*x):

f := (rs, z, xt, xs0, xs1) -> gga_exchange(pw91_f, rs, z, xs0, xs1):


text:
# Generated using bssn_4d_spher_matter.mw

restx:=expand(simplify((-2*b1(t, x)*alpha(t, x)*(diff(a1(t, x), t))-2*x*alpha(t, x)*(diff(a1(t, x), t))*(diff(b1(t, x), x))+2*x*(diff(diff(b1(t, x), t), x))*alpha(t, x)*a1(t, x)+2*(diff(b1(t, x), t))*alpha(t, x)*a1(t, x)+2*beta(t, x)*alpha(t, x)*b1(t, x)*(diff(a1(t, x), x))+b1(t, x)*x*alpha(t, x)*a1(t, x)*(diff(rpsi(t, x), x))*(diff(rpsi(t, x), t))+b1(t, x)*x*alpha(t, x)*a1(t, x)*(diff(ipsi(t, x), x))*(diff(ipsi(t, x), t))-b1(t, x)*x*alpha(t, x)*a1(t, x)*beta(t, x)*(diff(rpsi(t, x), x))^2-b1(t, x)*x*alpha(t, x)*a1(t, x)*beta(t, x)*(diff(ipsi(t, x), x))^2-2*beta(t, x)*alpha(t, x)*a1(t, x)*x*(diff(diff(b1(t, x), x), x))+2*beta(t, x)*alpha(t, x)*x*(diff(a1(t, x), x))*(diff(b1(t, x), x))-4*beta(t, x)*(diff(b1(t, x), x))*alpha(t, x)*a1(t, x)-2*x*(diff(b1(t, x), t))*a1(t, x)*(diff(alpha(t, x), x))+2*beta(t, x)*a1(t, x)*x*(diff(alpha(t, x), x))*(diff(b1(t, x), x))+2*beta(t, x)*a1(t, x)*b1(t, x)*(diff(alpha(t, x), x)))/(alpha(t, x)^3*b1(t, x)*x*a1(t, x))));

ire_rpsi_direct:=
-2*(diff(rpsi(t, x), x))*beta(t, x)*(diff(beta(t, x), x))/alpha(t, x)^2+(diff(rpsi(t, x), t))*beta(t, x)*(diff(a1(t, x), x))/(alpha(t, x)^2*a1(t, x))-(diff(rpsi(t, x), t))*beta(t, x)*(diff(alpha(t, x), x))/alpha(t, x)^3-(diff(rpsi(t, x), x))*beta(t, x)*(diff(alpha(t, x), t))/alpha(t, x)^3+(diff(rpsi(t, x), x))*beta(t, x)^2*(diff(alpha(t, x), x))/alpha(t, x)^3+2*(diff(rpsi(t, x), t))*beta(t, x)*(diff(b1(t, x), x))/(alpha(t, x)^2*b1(t, x))+2*(diff(rpsi(t, x), x))*beta(t, x)*(diff(b1(t, x), t))/(alpha(t, x)^2*b1(t, x))-2*(diff(rpsi(t, x), x))*beta(t, x)^2*(diff(b1(t, x), x))/(alpha(t, x)^2*b1(t, x))+(diff(rpsi(t, x), x))*beta(t, x)*(diff(a1(t, x), t))/(alpha(t, x)^2*a1(t, x))-(diff(rpsi(t, x), x))*beta(t, x)^2*(diff(a1(t, x), x))/(alpha(t, x)^2*a1(t, x))-2*(diff(rpsi(t, x), t))*(diff(b1(t, x), t))/(alpha(t, x)^2*b1(t, x))+2*(diff(rpsi(t, x), t))*beta(t, x)/(alpha(t, x)^2*x)-2*(diff(rpsi(t, x), x))*beta(t, x)^2/(alpha(t, x)^2*x)+(diff(rpsi(t, x), t))*(diff(alpha(t, x), t))/alpha(t, x)^3+2*beta(t, x)*(diff(diff(rpsi(t, x), t), x))/alpha(t, x)^2+(diff(rpsi(t, x), t))*(diff(beta(t, x), x))/alpha(t, x)^2-(diff(rpsi(t, x), t))*(diff(a1(t, x), t))/(alpha(t, x)^2*a1(t, x))-beta(t, x)^2*(diff(diff(rpsi(t, x), x), x))/alpha(t, x)^2+(diff(rpsi(t, x), x))*(diff(beta(t, x), t))/alpha(t, x)^2-(diff(diff(rpsi(t, x), t), t))/alpha(t, x)^2+2*(diff(rpsi(t, x), x))/(a1(t, x)^2*x)+(diff(diff(rpsi(t, x), x), x))/a1(t, x)^2-(diff(rpsi(t, x), x))*(diff(a1(t, x), x))/a1(t, x)^3+2*(diff(rpsi(t, x), x))*(diff(b1(t, x), x))/(b1(t, x)*a1(t, x)^2)-mass*rpsi(t, x)+(diff(rpsi(t, x), x))*(diff(alpha(t, x), x))/(alpha(t, x)*a1(t, x)^2);

resxx:=
(1/2)*(-alpha(t, x)^3*b1(t, x)^2*x^2*(diff(rpsi(t, x), x))^2-alpha(t, x)^3*b1(t, x)^2*x^2*(diff(ipsi(t, x), x))^2+alpha(t, x)*b1(t, x)^2*x^2*(diff(rpsi(t, x), x))^2*a1(t, x)^2*beta(t, x)^2+alpha(t, x)*b1(t, x)^2*x^2*(diff(ipsi(t, x), x))^2*a1(t, x)^2*beta(t, x)^2-2*b1(t, x)^2*a1(t, x)^2*alpha(t, x)*beta(t, x)^2-2*(diff(b1(t, x), t))^2*a1(t, x)^2*alpha(t, x)*x^2+4*b1(t, x)^2*x*alpha(t, x)^2*(diff(alpha(t, x), x))+4*(diff(b1(t, x), x))*alpha(t, x)^3*x*b1(t, x)+2*(diff(b1(t, x), x))^2*alpha(t, x)^3*x^2+2*alpha(t, x)^3*b1(t, x)^2-2*a1(t, x)^2*alpha(t, x)^3+alpha(t, x)^3*b1(t, x)^2*x^2*U(t, x)*a1(t, x)^2-alpha(t, x)*b1(t, x)^2*x^2*(diff(rpsi(t, x), t))^2*a1(t, x)^2-alpha(t, x)*b1(t, x)^2*x^2*(diff(ipsi(t, x), t))^2*a1(t, x)^2+4*a1(t, x)*beta(t, x)*b1(t, x)*x^2*alpha(t, x)*(diff(a1(t, x), t))*(diff(b1(t, x), x))-4*a1(t, x)*beta(t, x)^2*alpha(t, x)*b1(t, x)*x^2*(diff(a1(t, x), x))*(diff(b1(t, x), x))-4*a1(t, x)*beta(t, x)^2*alpha(t, x)*b1(t, x)^2*x*(diff(a1(t, x), x))+4*a1(t, x)*beta(t, x)*b1(t, x)^2*x*alpha(t, x)*(diff(a1(t, x), t))-4*b1(t, x)*x^2*a1(t, x)^2*alpha(t, x)*(diff(diff(b1(t, x), t), t))+4*b1(t, x)^2*x*a1(t, x)^2*alpha(t, x)*(diff(beta(t, x), t))-2*(diff(b1(t, x), x))^2*x^2*a1(t, x)^2*alpha(t, x)*beta(t, x)^2+4*b1(t, x)*x^2*(diff(alpha(t, x), t))*a1(t, x)^2*(diff(b1(t, x), t))-4*b1(t, x)^2*x*(diff(alpha(t, x), t))*a1(t, x)^2*beta(t, x)+4*b1(t, x)*x^2*alpha(t, x)^2*(diff(alpha(t, x), x))*(diff(b1(t, x), x))+4*b1(t, x)*x^2*a1(t, x)^2*alpha(t, x)*(diff(beta(t, x), t))*(diff(b1(t, x), x))+4*b1(t, x)*x^2*a1(t, x)^2*alpha(t, x)*beta(t, x)*(diff(diff(b1(t, x), t), x))+8*b1(t, x)*x*a1(t, x)^2*alpha(t, x)*beta(t, x)*(diff(b1(t, x), t))+4*(diff(b1(t, x), x))*x^2*a1(t, x)^2*alpha(t, x)*beta(t, x)*(diff(b1(t, x), t))-4*(diff(b1(t, x), x))*x*a1(t, x)^2*alpha(t, x)*beta(t, x)^2*b1(t, x)-4*alpha(t, x)*b1(t, x)^2*x*beta(t, x)*a1(t, x)^2*(diff(beta(t, x), x))-4*b1(t, x)*x^2*(diff(alpha(t, x), t))*a1(t, x)^2*beta(t, x)*(diff(b1(t, x), x))-4*alpha(t, x)*b1(t, x)*x^2*beta(t, x)*a1(t, x)^2*(diff(beta(t, x), x))*(diff(b1(t, x), x)))/(alpha(t, x)^3*a1(t, x)^2*b1(t, x)^2*x^2);

restt:=
(1/2)*(4*b1(t, x)^2*x*a1(t, x)^3*beta(t, x)^2*(diff(alpha(t, x), x))-2*(diff(b1(t, x), x))^2*x^2*a1(t, x)^3*alpha(t, x)*beta(t, x)^2+12*(diff(b1(t, x), x))*alpha(t, x)^3*a1(t, x)*x*b1(t, x)+4*alpha(t, x)^3*a1(t, x)*b1(t, x)*x^2*(diff(diff(b1(t, x), x), x))-4*alpha(t, x)^3*b1(t, x)*x^2*(diff(a1(t, x), x))*(diff(b1(t, x), x))-4*a1(t, x)^3*alpha(t, x)*b1(t, x)*x^2*beta(t, x)^2*(diff(diff(b1(t, x), x), x))+4*b1(t, x)*x^2*a1(t, x)^3*alpha(t, x)*beta(t, x)*(diff(diff(b1(t, x), t), x))+8*b1(t, x)*x*a1(t, x)^3*alpha(t, x)*beta(t, x)*(diff(b1(t, x), t))+4*(diff(b1(t, x), x))*x^2*a1(t, x)^3*alpha(t, x)*beta(t, x)*(diff(b1(t, x), t))-12*(diff(b1(t, x), x))*x*a1(t, x)^3*alpha(t, x)*beta(t, x)^2*b1(t, x)+4*a1(t, x)^3*(diff(b1(t, x), t))*alpha(t, x)*b1(t, x)*x^2*(diff(beta(t, x), x))-4*a1(t, x)^2*(diff(b1(t, x), t))*alpha(t, x)*b1(t, x)*x^2*(diff(a1(t, x), t))-4*alpha(t, x)*b1(t, x)^2*x*beta(t, x)*a1(t, x)^3*(diff(beta(t, x), x))-4*b1(t, x)*x^2*a1(t, x)^3*(diff(b1(t, x), t))*beta(t, x)*(diff(alpha(t, x), x))+4*b1(t, x)*x^2*a1(t, x)^3*beta(t, x)^2*(diff(b1(t, x), x))*(diff(alpha(t, x), x))+4*a1(t, x)^2*(diff(b1(t, x), t))*alpha(t, x)*b1(t, x)*x^2*beta(t, x)*(diff(a1(t, x), x))-4*alpha(t, x)*b1(t, x)*x^2*beta(t, x)*a1(t, x)^3*(diff(beta(t, x), x))*(diff(b1(t, x), x))-2*b1(t, x)^2*a1(t, x)^3*alpha(t, x)*beta(t, x)^2-2*(diff(b1(t, x), t))^2*a1(t, x)^3*alpha(t, x)*x^2+2*(diff(b1(t, x), x))^2*alpha(t, x)^3*a1(t, x)*x^2-4*alpha(t, x)^3*b1(t, x)^2*x*(diff(a1(t, x), x))+2*alpha(t, x)^3*a1(t, x)*b1(t, x)^2-2*a1(t, x)^3*alpha(t, x)^3+b1(t, x)^2*x^2*alpha(t, x)*a1(t, x)^3*(diff(rpsi(t, x), t))^2+b1(t, x)^2*x^2*alpha(t, x)*a1(t, x)^3*(diff(ipsi(t, x), t))^2+b1(t, x)^2*x^2*alpha(t, x)^3*a1(t, x)^3*U(t, x)+b1(t, x)^2*x^2*alpha(t, x)^3*a1(t, x)*(diff(rpsi(t, x), x))^2+b1(t, x)^2*x^2*alpha(t, x)^3*a1(t, x)*(diff(ipsi(t, x), x))^2-b1(t, x)^2*x^2*alpha(t, x)*a1(t, x)^3*(diff(rpsi(t, x), x))^2*beta(t, x)^2-b1(t, x)^2*x^2*alpha(t, x)*a1(t, x)^3*(diff(ipsi(t, x), x))^2*beta(t, x)^2)/(alpha(t, x)^3*b1(t, x)^2*x^2*a1(t, x)^3);

resthth:=
(1/2)*(-2*a1(t, x)^2*alpha(t, x)*b1(t, x)*x*(diff(diff(a1(t, x), t), t))+2*alpha(t, x)^2*a1(t, x)*b1(t, x)*x*(diff(diff(alpha(t, x), x), x))-2*alpha(t, x)*a1(t, x)^3*b1(t, x)*x*(diff(beta(t, x), x))^2-2*b1(t, x)*x*(diff(a1(t, x), x))*alpha(t, x)^2*(diff(alpha(t, x), x))+2*(diff(alpha(t, x), x))*a1(t, x)^2*b1(t, x)*x*beta(t, x)^2*(diff(a1(t, x), x))+2*(diff(alpha(t, x), x))*a1(t, x)^3*b1(t, x)*x*beta(t, x)*(diff(beta(t, x), x))-2*(diff(alpha(t, x), x))*a1(t, x)^2*b1(t, x)*x*beta(t, x)*(diff(a1(t, x), t))+2*a1(t, x)^2*(diff(b1(t, x), t))*alpha(t, x)*x*beta(t, x)*(diff(a1(t, x), x))-2*alpha(t, x)*x*beta(t, x)^2*a1(t, x)^2*(diff(a1(t, x), x))*(diff(b1(t, x), x))-4*alpha(t, x)*x*beta(t, x)*a1(t, x)^3*(diff(beta(t, x), x))*(diff(b1(t, x), x))+2*alpha(t, x)*x*beta(t, x)*a1(t, x)^2*(diff(a1(t, x), t))*(diff(b1(t, x), x))-2*a1(t, x)^2*(diff(alpha(t, x), t))*b1(t, x)*x*beta(t, x)*(diff(a1(t, x), x))+4*(diff(a1(t, x), t))*alpha(t, x)*a1(t, x)^2*b1(t, x)*x*(diff(beta(t, x), x))+2*a1(t, x)^2*alpha(t, x)*b1(t, x)*x*(diff(beta(t, x), t))*(diff(a1(t, x), x))+4*a1(t, x)^2*alpha(t, x)*b1(t, x)*x*beta(t, x)*(diff(diff(a1(t, x), t), x))-2*alpha(t, x)*a1(t, x)^2*b1(t, x)*x*beta(t, x)^2*(diff(diff(a1(t, x), x), x))-2*alpha(t, x)*a1(t, x)^3*b1(t, x)*x*beta(t, x)*(diff(diff(beta(t, x), x), x))-6*alpha(t, x)*a1(t, x)^2*b1(t, x)*x*beta(t, x)*(diff(a1(t, x), x))*(diff(beta(t, x), x))+2*b1(t, x)*a1(t, x)^3*beta(t, x)^2*(diff(alpha(t, x), x))+4*a1(t, x)^3*alpha(t, x)*beta(t, x)*(diff(b1(t, x), t))-4*(diff(b1(t, x), x))*a1(t, x)^3*alpha(t, x)*beta(t, x)^2-2*x*a1(t, x)^3*alpha(t, x)*(diff(diff(b1(t, x), t), t))+2*b1(t, x)*a1(t, x)^3*alpha(t, x)*(diff(beta(t, x), t))+2*a1(t, x)*b1(t, x)*alpha(t, x)^2*(diff(alpha(t, x), x))+2*x*(diff(alpha(t, x), t))*a1(t, x)^3*(diff(b1(t, x), t))-2*b1(t, x)*(diff(alpha(t, x), t))*a1(t, x)^3*beta(t, x)+2*alpha(t, x)^3*a1(t, x)*x*(diff(diff(b1(t, x), x), x))-2*alpha(t, x)^3*x*(diff(a1(t, x), x))*(diff(b1(t, x), x))+b1(t, x)*x*a1(t, x)^3*alpha(t, x)^3*U(t, x)-b1(t, x)*x*a1(t, x)^3*alpha(t, x)*(diff(rpsi(t, x), t))^2-b1(t, x)*x*a1(t, x)^3*alpha(t, x)*(diff(ipsi(t, x), t))^2+b1(t, x)*x*a1(t, x)*alpha(t, x)^3*(diff(rpsi(t, x), x))^2+b1(t, x)*x*a1(t, x)*alpha(t, x)^3*(diff(ipsi(t, x), x))^2-2*alpha(t, x)^3*b1(t, x)*(diff(a1(t, x), x))+4*(diff(b1(t, x), x))*alpha(t, x)^3*a1(t, x)+2*a1(t, x)^3*(diff(b1(t, x), t))*alpha(t, x)*x*(diff(beta(t, x), x))-2*a1(t, x)^2*(diff(b1(t, x), t))*alpha(t, x)*x*(diff(a1(t, x), t))-2*alpha(t, x)*b1(t, x)*beta(t, x)^2*a1(t, x)^2*(diff(a1(t, x), x))-4*alpha(t, x)*b1(t, x)*beta(t, x)*a1(t, x)^3*(diff(beta(t, x), x))+2*alpha(t, x)*b1(t, x)*beta(t, x)*a1(t, x)^2*(diff(a1(t, x), t))-2*a1(t, x)^3*(diff(alpha(t, x), t))*b1(t, x)*x*(diff(beta(t, x), x))+2*a1(t, x)^2*(diff(alpha(t, x), t))*b1(t, x)*x*(diff(a1(t, x), t))+2*a1(t, x)^3*alpha(t, x)*b1(t, x)*x*(diff(diff(beta(t, x), t), x))-2*a1(t, x)^3*alpha(t, x)*x*beta(t, x)^2*(diff(diff(b1(t, x), x), x))+2*x*a1(t, x)^3*alpha(t, x)*(diff(beta(t, x), t))*(diff(b1(t, x), x))+4*x*a1(t, x)^3*alpha(t, x)*beta(t, x)*(diff(diff(b1(t, x), t), x))+2*a1(t, x)*x*alpha(t, x)^2*(diff(alpha(t, x), x))*(diff(b1(t, x), x))-2*x*(diff(alpha(t, x), t))*a1(t, x)^3*beta(t, x)*(diff(b1(t, x), x))-2*x*a1(t, x)^3*(diff(b1(t, x), t))*beta(t, x)*(diff(alpha(t, x), x))+2*x*a1(t, x)^3*beta(t, x)^2*(diff(b1(t, x), x))*(diff(alpha(t, x), x))-b1(t, x)*x*a1(t, x)^3*alpha(t, x)*(diff(rpsi(t, x), x))^2*beta(t, x)^2-b1(t, x)*x*a1(t, x)^3*alpha(t, x)*(diff(ipsi(t, x), x))^2*beta(t, x)^2+2*b1(t, x)*x*a1(t, x)^3*alpha(t, x)*beta(t, x)*(diff(rpsi(t, x), t))*(diff(rpsi(t, x), x))+2*b1(t, x)*x*a1(t, x)^3*alpha(t, x)*beta(t, x)*(diff(ipsi(t, x), t))*(diff(ipsi(t, x), x)))/(b1(t, x)*x*alpha(t, x)^3*a1(t, x)^3);


text:
is_left_unit := (A::set) -> proc(m,e) local a;
 `and`(seq(evalb(m(e,a)=a),a in A)); end:

is_right_unit := (A::set) -> proc(m,e) local a;
 `and`(seq(evalb(m(a,e)=a),a in A)); end:

is_left_absorbing := (A::set) -> proc(m,z) local a;
 `and`(seq(evalb(m(z,a)=z),a in A)); end:

is_right_absorbing := (A::set) -> proc(m,z) local a;
 `and`(seq(evalb(m(a,z)=z),a in A)); end:

is_idempotent := (A::set) -> proc(m) local a;
 `and`(seq(evalb(m(a,a)=a),a in A)); end:

is_commutative := (A::set) -> proc(m) local a,b;
 `and`(seq(seq(evalb(m(a,b) = m(b,a)),b in A),a in A)); end:

is_associative := (A::set) -> proc(m) local a,b,c;
 `and`(seq(seq(seq(evalb(m(a,m(b,c))=m(m(a,b),c)),c in A),b in A),a in A)); end:

is_left_distributive := (A::set) -> proc(m,p) local a,b,c;
 `and`(seq(seq(seq(evalb(m(a,p(b,c))=p(m(a,b),m(a,c))),c in A),b in A),a in A)); end:

is_right_distributive := (A::set) -> proc(m,p) local a,b,c;
 `and`(seq(seq(seq(evalb(m(p(a,b),c)=p(m(a,c),m(b,c))),c in A),b in A),a in A)); end:

is_monoid := (A::set) -> proc(m,e)
 is_left_unit(A)(m,e) and is_right_unit(A)(m,e) and is_associative(A)(m); end:

is_semiring := (A::set) -> (m,p,e,z) ->
 is_monoid(A)(m,e) and
 is_monoid(A)(p,z) and
 is_left_absorbing(A)(m,z) and
 is_right_absorbing(A)(m,z) and
 is_left_distributive(A)(m,p) and
 is_right_distributive(A)(m,p);


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_c *)

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

malpha := 2.804:
mgamma := 0.8098:

XX := s -> 1/(1 + malpha*s^2):
ff := s -> XX(s) + mgamma*(1 - XX(s)):

f := (rs, z, xt, xs0, xs1) -> f_pw(rs, z)*(
  + (1 + z)/2 * ff(X2S*xs0)
  + (1 - z)/2 * ff(X2S*xs1)
):


text:
#!/bin/bash maple
# Use `maple -q symengine_bench.mpl -D n=15` to run

e := sin(cos(x+1)):
st := time[real]():
f := series(e, x=0,n):
1000*(time[real]()-st);

done

max_issues_repo_name: HQSquantumsimulations/symengine

max_forks_repo_name: HQSquantumsimulations/symengine


text:
program assignchar;	{sample23}
var x,y : char;;
begin
 x := 'X';
 y := 'Y'
end.

max_issues_repo_name: shouth/LanguageProcessing


text:
`depends/Bucket` := proc(mr, r::(name=range), x, $)
  depends(rhs(r), x) or depends(mr, x minus {lhs(r)})
end proc:

`depends/Index` := proc(n, o::name, e, mr, x, $)
  depends({n,e}, x) or depends(mr, x minus {o})
end proc:

# note that v _can_ occur in m1.
`depends/Let` := proc(m1, v::name, m2, x, $)
  depends(m1, x) or depends(m2, x minus {v})
end proc:

`eval/Bucket` := proc(e::anyfunc(anything,name=range), eqs, $)
  local bvar, body;
  bvar, body := BindingTools:-generic_evalat(op([2,1],e), op(1,e), eqs);
  eval(op(0,e), eqs)(body, bvar = eval(op([2,2],e), eqs))
end proc:

`eval/Index` := proc(e::anyfunc(anything,name,anything,anything), eqs, $)
  local o, mr;
  o, mr := BindingTools:-generic_evalat(op(2,e), op(4,e), eqs);
  eval(op(0,e), eqs)(eval(op(1,e), eqs), o, eval(op(3,e), eqs), mr)
end proc:

`eval/Let` := proc(e, eqs, $)
  local m1, v, m2;
  m1, v, m2 := op(e);
  eval(op(0,e), eqs)(eval(m1, eqs), BindingTools:-generic_evalat(v, m2, eqs))
end proc:

Summary := module ()
  option package;
  export RoundTrip, Summarize, SummarizeKB, bucket, summarize;
  global Bucket, Fanout, Index, Split, Nop, Add, Let,
         SumIE, ProductIE, BucketIE;
  uses Hakaru, KB;

  RoundTrip := proc(e)
    local result;
    interface(screenwidth=9999, prettyprint=0, warnlevel=0,
      showassumed=0,quiet=true);
    kernelopts(assertlevel=0);
    result := eval(ToInert(Summarize(_passed)), _Inert_ATTRIBUTE=NULL);
    lprint(result)
  end proc;

  Summarize := proc(e, {ctx :: list := []}, $)
    local ie;
    ie := table('[sum    =SumIE    , Sum    =SumIE    ,
                  product=ProductIE, Product=ProductIE,
                  bucket =BucketIE , Bucket =BucketIE ]');
    subsindets(SummarizeKB(e, foldr(assert, empty, op(ctx))),
               And(specfunc({indices(ie, 'nolist')}),
                   anyfunc(anything, anything=range)),
               e -> subsop(0 = ie[op(0,e)],
                           applyop(`+`, [2,2,2], e, 1)))
  end proc;

  SummarizeKB := proc(e, kb :: t_kb, $)
    local patterns, x, kb1, e1, rng, summary, mr, f;
    if not hasfun(e, '{Sum,sum}', 'piecewise') then
      e;
    elif e :: '{known_continuous, known_discrete,
                specfunc({Msum, Weight, Fanout, Split, Nop, Add,
                          exp, log, GAMMA, Beta, idx, And, Or, Not}),
                `+`, `*`, `^`, `..`, boolean, indexed, list}' then
      map(procname, _passed);
    elif e :: 'Context(anything, anything)' then
      kb1 := assert(op(1,e), kb);
      applyop(procname, 2, e, kb1);
    elif e :: 'specfunc(piecewise)' then
      kb_piecewise(e, kb, ((lhs, kb) -> lhs), SummarizeKB);
    elif e :: 'lam(name, anything, anything)' then
      patterns := htype_patterns(op(2,e));
      if patterns :: Branches(Branch(PVar(name),anything)) then
        # Eta-expand the function type
        x := op(1,e);
        x, kb1 := genType(x, op(2,e), kb);
        e1 := eval(op(3,e), op(1,e)=x);
        lam(x, op(2,e), SummarizeKB(e1, kb1))
      else
        # Eta-expand the function type and the sum-of-product argument-type
        x := op(1,e);
        if depends([e,kb], x) then x := gensym(x) end if;
        e1 := eval(op(3,e), op(1,e)=x);
        lam(x, op(2,e), 'case'(x,
          map(proc(branch)
                local eSubst, pSubst, p1, binds, y, kb1, i, pSubst1;
                eSubst, pSubst := pattern_match([x,e], x, op(1,branch));
                p1 := subs(pSubst, op(1,branch));
                binds := [pattern_binds(p1)];
                kb1 := kb;
                pSubst1 := table();
                for i from 1 to nops(binds) do
                  y, kb1 := genType(op(i,binds), op([2,i],branch), kb1);
                  pSubst1[op(i,binds)] := y;
                end do;
                pSubst1 := op(op(pSubst1));
                Branch(subs(pSubst1, p1),
                       SummarizeKB(eval(eval(e1,eSubst),pSubst1), kb1))
              end proc,
              patterns)))
      end if
    elif e :: 'ary(anything, name, anything)' then
      x, kb1 := genType(op(2,e), HInt(closed_bounds(0..op(1,e)-1)), kb);
      e1 := SummarizeKB(eval(op(3,e), op(2,e)=x), kb1);
      subsop(2=x, 3=e1, e);
    elif e :: 'And(specfunc({sum, Sum, product, Product}),
                   anyfunc(anything, name=range))' then
      rng := op([2,2],e);
      x, kb1 := genType(op([2,1],e), HInt(closed_bounds(rng)), kb);
      rng := map(SummarizeKB, rng, kb);
      e1 := eval(op(1,e), op([2,1],e)=x);
      if op(0, e) in '{sum, Sum}' and has(e1, 'piecewise') then
        mr, f := summarize(e1, kb, x, summary);
        if hasfun(mr, '{Fanout, Index}') then
          mr := SummarizeKB(mr, kb1);
          return Let(Bucket(mr, x=rng),
                     summary,
                     NewSLO:-simplify_factor_assuming(f, kb));
        end if;
      end if;
      e1 := SummarizeKB(e1, kb1);
      subsop(2=(x=rng), 1=e1, e);
    else
      error "SummarizeKB doesn't know how to handle %1", e;
    end if;
  end proc;

  bucket := proc(mr, r::(name=range), cond::list:=[], $)
    if mr :: 'Fanout(anything,anything)' then
      Pair(op(map(procname, _passed)))
    elif mr :: 'Split(anything,anything,anything)' then
      Pair(bucket(op(2,mr), r, [    op(1,mr) ,op(cond)]),
           bucket(op(3,mr), r, [Not(op(1,mr)),op(cond)]))
    elif mr :: 'Index(anything,name,anything,anything)' then
      ary(op(1,mr), op(2,mr),
          bucket(op(4,mr), r, [op(2,mr)=op(3,mr),op(cond)]));
    elif mr :: 'Nop()' then
      _Unit
    elif mr :: 'Add(anything)' then
      sum(piecewise_And(cond, op(1,mr), 0), r)
    else
      'procname(_passed)'
    end if;
  end proc;

  summarize := proc(ee,
                    kb :: t_kb,
                    i :: {name,list(name),set(name)},
                    summary, $)
    local e, r, s,
          e1, mr1, f1, e2, mr2, f2,
          variables, var_outerness, outermost, thunk, consider,
          o, t, lo, hi, b, a;

    # Try to ensure termination
    e := simplify_assuming(ee, kb);
    if length(e) >= length(ee) then e := ee end if;

    # Nop rule
    if Testzero(e) then return Nop(), 0 end if;

    r := indets(e, '{relation, logical, specfunc({And,Or})}');
    s, r := selectremove(depends, r, i);
    # Fanout rule
    r := sort(convert(r, 'list'), 'length');
    while nops(r) > 0 do
      e1 := eval(e, r[-1]=true);  if e = e1 then r := r[1..-2]; next end if;
      e2 := eval(e, r[-1]=false); if e = e2 then r := r[1..-2]; next end if;
      mr1, f1 := summarize(e1, kb, i, 'fst(summary)');
      mr2, f2 := summarize(e2, kb, i, 'snd(summary)');
      return Fanout(mr1, mr2), 'piecewise'(r[-1], f1, f2);
    end do;

    variables     := kb_to_variables(kb);
    var_outerness := table(variables =~ [seq(1..nops(variables))]);
    outermost     := -infinity;
    consider      := proc(expr, th, $)
      local outerness;
      outerness := map((n -> var_outerness[n]), indets(expr, 'name'));
      outerness := min(select(type, outerness, 'integer'));
      if outermost < outerness then
        outermost := outerness;
        thunk     := th;
      end if;
    end proc;
    if has(variables, i) then
      error "The loop to be summarized shadows the variable %1", i;
    end if;
    for r in s do
      e1 := eval(e, r=true);  if e = e1 then next end if;
      e2 := eval(e, r=false); if e = e2 then next end if;
      # Index rule
      if r :: `=` and Testzero(e2) then
        for o in indets(r, 'name') minus convert(i, 'set') do
          if not (var_outerness[o] :: integer) then next end if;
          t := getType(kb, o);
          if not (t :: specfunc(HInt)) then next end if;
          lo := op(select(type, t, 'Bound(identical(`>=`), anything)'));
          if lo = NULL then next else lo := op(2,lo) end if;
          hi := op(select(type, t, 'Bound(identical(`<=`), anything)'));
          if hi = NULL then next else hi := op(2,hi) end if;
          if not ispoly(`-`(op(r)), 'linear', o, 'b', 'a') then next end if;
          b := Normalizer(-b/a);
          consider(b, ((e1, o, lo, hi, b) -> proc($)
            local mr, f;
            mr, f := summarize(e1, kb, i, 'idx'(summary, o-lo));
            Index(hi-lo+1, o, b, mr),
              'piecewise'(And(o::integer, lo<=o, o<=hi), f, 0);
          end proc)(e1, o, lo, hi, b));
        end do;
      end if;
      # Split rule
      consider(r, ((e1, e2, r) -> proc($)
        local mr1, f1, mr2, f2;
        mr1, f1 := summarize(e1, kb, i, 'fst(summary)');
        mr2, f2 := summarize(e2, kb, i, 'snd(summary)');
        Split(r, mr1, mr2), f1 + f2;
      end proc)(e1, e2, r));
    end do;
    if thunk :: procedure then
      return thunk();
    else
      # Add rule
      return Add(e), summary;
    end if;
  end proc;
end module; # Summary

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


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 *)

mkappa := 0.4604:
mmu    := 0.354546875:

f0 := s -> 1 + mkappa*(1 - mkappa/(mkappa + mmu*s^2)):
f  := x -> f0(X2S_2D*x):


text:
#@ Not autoload

with(GroupTheory):
with(LinearAlgebra):

S4 := SymmetricGroup(4);

x2 := [[1,2]];
y2 := [[2,3]];
z2 := [[3,4]];
u2 := [[1,2],[3,4]];
v2 := [[1,3],[2,4]];
x3 := [[1,2,3]];
x4 := [[1,3,2,4]];

H[ 1] := Subgroup({},S4);
H[ 2] := Subgroup({x2},S4);
K[ 2] := Subgroup({u2},S4);
H[ 3] := Subgroup({x3},S4);
H[ 4] := Subgroup({x4},S4);
K[ 4] := Subgroup({x2,z2},S4);
L[ 4] := Subgroup({u2,v2},S4);
H[ 6] := Subgroup({x2,x3},S4);
H[ 8] := Subgroup({x2,x4},S4);
H[12] := Subgroup({u2,x3},S4);
H[24] := S4;

S4_subgroups := [H[1],H[2],K[2],H[3],H[4],K[4],L[4],H[6],H[8],H[12],H[24]]:

# Burnside matrix of numbers |Map_G(G/H_i,G/H_j)|
n := nops(S4_subgroups);
M := Matrix(n,n):
for i from 1 to n do
 for j from 1 to n do
  Hi := S4_subgroups[i];
  Zi := {op(Generators(Hi))};
  Hj := S4_subgroups[j];
  C := map(Representative,LeftCosets(Hj,S4));
  M[i,j] := nops(select(c -> `and`(op(map(x -> SubgroupMembership(c^(-1) . x . c,Hj),Zi))),C));
 od:
od:

unassign('C'):

# Skeleton of the poset of conjugacy classes of subgroups
N := {seq(i,i=1..n)};
R := select(ij -> M[ij[1],ij[2]] <> 0,{seq(seq([i,j],j=1..n),i=1..n)});

# Conjugacy classes

C[1] := ConjugacyClass([],S4);
C[2] := ConjugacyClass(x2,S4);
C[3] := ConjugacyClass(u2,S4);
C[4] := ConjugacyClass(x3,S4);
C[5] := ConjugacyClass(x4,S4);

chars := [
 [ 1, 1, 1, 1, 1],
 [ 1,-1, 1, 1,-1],
 [ 2, 0, 2,-1, 0],
 [ 3, 1,-1, 0,-1],
 [ 3,-1,-1, 0, 1]
];

char_names := [1,epsilon,sigma,rho,epsilon*rho];

char_rels := [
 epsilon^2 - 1,
 epsilon * sigma - sigma,
 sigma^2 - 1 - epsilon - sigma,
 sigma * rho - rho - epsilon * rho,
 rho^2 - 1 - sigma - rho - epsilon * rho
];

char_psi_rels := [
 psi(2,epsilon) = 1,
 psi(2,sigma) = 1 - epsilon + sigma,
 psi(3,sigma) = 1 + epsilon,
 psi(2,rho) = 1 + sigma + rho - epsilon * rho,
 psi(3,rho) = 1 + epsilon - sigma + rho
];

char_lambda_rels := [
 lambda(2,sigma) = epsilon,
 lambda(2,rho) = epsilon * rho,
 lambda(3,rho) = epsilon
];


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: mgga_exc *)

$define mgga_x_revtpss_params
$include "mgga_x_tpss.mpl"

reg_c := 3:
reg_d := 1.475:

reg_f_a :=  a -> (1 - a)^3/(1 + (reg_d*a)^2)^(3/2):

(* Eq. (12). Note that alpha = 0 => t = x^2/8 *)
reg_f := (x, u, t) ->
  tpss_f(x, u, t) + reg_f_a(tpss_alpha(x, t))*exp(-reg_c*X2S^2*x^2)*(tpss_f(x, u, x^2/8) - tpss_f(x, u, t)):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) ->
  mgga_exchange(reg_f, rs, z, xs0, xs1, u0, u1, t0, t1):


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 *)
(* prefix:
  gga_x_vmt_params *params;

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

f0_vmt := s -> 1 + params_a_mu*s^2*exp(-params_a_alpha*s^2)/(1 + params_a_mu*s^2):

f_vmt  := x -> f0_vmt(X2S*x):
f      := x -> f_vmt(x):


text:
   {   Transp.mpl   }

TITLE
    TransportationProblem;

INDEX
    supply = (Wh1,Wh2,Wh3,Wh4,Wh5);
    dest   = (A1,A2,B1,B2,B3,C1,C2,C3);

DATA
    MaxCapacity[supply]  := (200,300,250,450,350);
    Required[dest]       := (188,214,137,199,225,200,177,185);

    ShippingCost[supply,dest] := (13, 22,  9, 14, 10, 32, 35, 36,
                                  12, 10, 11, 13,  7, 11, 16, 15,
                                  24, 22, 17, 11,  9,  8, 12,  8,
                                  43, 40, 25, 17, 21, 24, 27, 22,
                                  15,  9, 41, 38, 35,  8,  5, 13);

DECISION VARIABLES
    VolumeShipped[supply,dest]  ->  ""

MODEL

    MIN  TotalCost = SUM(supply,dest: ShippingCost * VolumeShipped);

SUBJECT TO

    Capacity[supply] : SUM(dest: VolumeShipped)    <  MaxCapacity ;
    Demand[dest]     : SUM(supply: VolumeShipped)  >  Required ;

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models