print(`Sophia for Maple R6 - 13 june 2000, based on version 2 December 1998`); ############################################################################## # DYNAMIC TOOLS FOR MAPLE # Sophia # By # Martin Lesser and Anders Lennartsson # Royal Institute of Technology # Stockholm, Sweden # © 1992-1998 # M.B. Lesser and Anders Lennartsson ############################################################################## # DATA TYPES and GLOBAL VARIABLES ############################################################################## # Rotation Matrices # Evectors # Edyads # Kvectors # KMvectors # SKvectors # Partitioned vectors ############################################################################## # Frame names as they are created e.g. RAB ############################################################################## # Sets for making various substitutions ############################################################################## # toTimeFunction # toTimeExpression # SetOfFramesUsed # kde ############################################################################## # Loading of needed packages ############################################################################## with(linalg): ############################################################################## # Rotation Matrices, Creation Modification and Accounting ############################################################################## # A frame relation is a set containing the names of two frames # all these sets are stored in the global set SetOfFramesUsed # this is used to record all direction cosine matrices known # to the system. ############################################################################## declareFrameRelation:= proc(Frz1,Frz2) global SetOfFramesUsed; if not assigned(SetOfFramesUsed) then SetOfFramesUsed:={} fi; SetOfFramesUsed := {op(SetOfFramesUsed),{Frz1,Frz2}}: end: ############################################################################## # Sets up a standard Maple 3by3 array, also records this # fact in the set SetOfFramesUsed. These direction cosine # matrices are defined with names starting with R followed by # the frame names. Two forms of input are dealt with. If the # call contains 11 arguments the first two are frame names and # the last 9 are direction cosines. If only three args are use # they are all interpreted as frame names, the third being # an intermediate frame, say C. RAB and RBA are then computed # from the known RAC, RCB. ############################################################################## # defines the dircos matrix for a simple rotation consisting of # a rotation in the positive sense about an axis. The positive # sense of a rotation follows the usual cross product rule ############################################################################## SimpRot := proc(A,B,axis,ang) SophiaDeclareFrameRelation(B,[A,axis,ang]); end: ############################################################################## # Up to four angles (five frames) are related by simple rotations # about coordinate axis of the proceeding case. This is input # in the form of a list of lists. The component lists have the # same form as used in the simprot procedure. ############################################################################## chainSimpRot := proc() local all,each; if type([args],list(STRotationList)) then all:=[args]; for each in all do SophiaDeclareFrameRelation(each[2],[each[1],each[3],each[4]]); od; elif type([args],[list(STRotationList)]) then all:=args; for each in all do SophiaDeclareFrameRelation(each[2],[each[1],each[3],each[4]]); od; else ERROR(`Invalid argument(s)!`); fi; end: ############################################################################## # Just a rehash of simpRot with a better convention # for the order of frames in the arguments # as in many other cases older forms are kept to allow older # files to work ############################################################################## rotFromTo := proc(A,B,axis,ang) SophiaDeclareFrameRelation(B,[A,axis,ang]); end: ############################################################################## # A monodic operator form of the simp rot procedure ############################################################################## `&rot` := proc(rotList) rotFromTo(op(rotList)) end: ############################################################################## # EVECTORS AND EDYADS CONSTRUCTION AND SELECTION ############################################################################## # The Evector is just a list of an array and a symbol or # frame name. It is a way of keeping track of in which frame # a vector is represented. This is an Evector constructor ############################################################################## Evector := proc(v1,v2,v3,frameName) RETURN(SophiaConstructEvector([v1,v2,v3],frameName)); end: ############################################################################## # The same for Dyads or second rank tensors ############################################################################## Edyad:=proc(d11,d12,d13,d21,d22,d23,d31,d32,d33,frameName) if type([args],[algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,STuaname]) then RETURN([array(1..3,1..3,[[d11,d12,d13],[d21,d22,d23],[d31,d32,d33]]),frameName]); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # an dyadic operator form of Edyad constructor from a frame # name and an array rather than a list of components ############################################################################## `&Dyad` := proc(frm,ary) [ary,frm]; end: ############################################################################## # operator form of constructor. left arg is frame name and # the right arg is a list of components in that frame ############################################################################## `&ev` := proc(Fr,ve) Evector(ve[1],ve[2],ve[3],Fr) end: ############################################################################## # A way of constructing unit vectors in base vector directions # thus A&>1 is the unit vector a1 for frame A. ############################################################################## `&>` := proc(frm,nz) if nz=1 then RETURN(Evector(1,0,0,frm)); elif nz=2 then RETURN(Evector(0,1,0,frm)); elif nz=3 then RETURN(Evector(0,0,1,frm)); fi; end: ############################################################################## # selectors for the vector part (array) or frame part (name) # of an Evector or Edyad (gives 3 by 3 array) ############################################################################## `&vPart` := proc(ev) RETURN(ev[1]); end: `&fPart` := proc(ev) RETURN(ev[2]); end: ############################################################################## # selects component of an Evector in whatever frame it is # represented, i.e. the frame (&fPart ev). ############################################################################## Ec:=proc(evec,axi) RETURN(evec[1][axi]); end: ############################################################################## # A dyadic operator that gives the component ############################################################################## `&c` := proc(ev,ax) RETURN(Ec(ev,ax)); end: ############################################################################## # TRANSFORMATIONS OF EVECTORS AND EDYADS ############################################################################## # Applies the maple simplifier to components of Evectors and # Edyads. Beware of its primative construction not using # buffered construction ############################################################################## Esimplify := proc(evec) RETURN([map(simplify,&vPart evec),&fPart evec]); end: ############################################################################## # transforms a dyad from its present frame to a new frame ############################################################################## expressDyad := proc(edyad,newF) local oldFrame,RNewOld,ROldNew; oldFrame := (&fPart edyad); if oldFrame=newF then RETURN(edyad) fi; RNewOld := SophiaSelectRmx(newF,oldFrame); ROldNew := SophiaSelectRmx(oldFrame,newF); RETURN([evalm(RNewOld &* (&vPart edyad) &* ROldNew),newF]); end: ############################################################################## # transforms components of a dyad or a vector between frames ############################################################################## #type(DI[1],'matrix'(algebraic,square)); ############################################################################## express := proc(evec,newF) local oldFrame,RNewOld,ROldNew; oldFrame := (&fPart evec); if oldFrame=newF then RETURN(evec) fi; RNewOld := SophiaSelectRmx(newF,oldFrame); # added by Lennartsson to make express work on dyads also if type(evec[1],'matrix'(algebraic,square)) then ROldNew := SophiaSelectRmx(oldFrame,newF); RETURN([evalm(RNewOld &* (&vPart evec) &* ROldNew),newF]); else RETURN([evalm(RNewOld &* (&vPart evec)),newF]); fi; end: ############################################################################## # combines transformation with simplification operation ############################################################################## sexpress:=proc(Evec,fram) RETURN(Esimplify(express(Evec,fram))); end: ############################################################################## # a dyadic operator form of transformation ############################################################################## `&to` := proc(Fr,ev) RETURN(sexpress(ev,Fr)); end: ############################################################################## # monodic operator for simplification of Evectors and Edyads ############################################################################## `&simp` := proc(ev) RETURN(Esimplify(ev)); end: ############################################################################## # provides a form with trig functions as sums of angles ############################################################################## EsimplifyTrig := proc(evec) RETURN([map(combine,map(simplify, &vPart evec),trig), &fPart evec]); end: ############################################################################## # ALGEBRA OF EVECTORS AND DYADS ############################################################################## ############################################################################## # dyad dotted into Evector produces an Evector, includes # simplification. Needs modification as it does not use # a constructor as buffer to Edyad form ############################################################################## EdyadDot:= proc(Edyd,Evec) local ev1; if ( &fPart Edyd)= ( &fPart Evec) then RETURN([multiply( &vPart Edyd, &vPart Evec), &fPart Edyd]) fi; ev1:= sexpress(Evec, &fPart Edyd); Esimplify([multiply( &vPart Edyd, &vPart ev1), &fPart ev1]); end: ############################################################################## # ############################################################################## Edot:= proc(ej1,ej2) local ej3; if type([args],[STEvector,STEvector]) then ej3:=express(ej2,ej1[2]); RETURN(dotprod(ej1[1],ej3[1],orthogonal)); elif type([args],[STEdyad,STEvector]) then ej3:=express(ej2,ej1[2]); RETURN([multiply(ej1[1],ej3[1]),ej1[2]]); elif type([args],[STEvector,STEdyad]) then ej3:=express(ej2,ej1[2]); RETURN([multiply(transpose(ej3[1]),ej1[1]),ej1[2]]); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # dyadic operator form of Evector dot product ############################################################################## `&o` := proc(ev1,ev2) RETURN(Edot(ev1,ev2)); end: ############################################################################## # adds form of two evectors multiplied by scalars ############################################################################## EaddMixed := proc(s1,evec1,s2,evec2,ansFrame) local ev1,ev2; ev1 := express(evec1,ansFrame); ev2 := express(evec2,ansFrame); Eadd(s1,ev1,s2,ev2) end: ############################################################################## # adds up a list of Evectors ############################################################################## EaddList := proc(EvecList,ansFrame) local Eresult, Ev; if type(EvecList[1][1],matrix) then Eresult := Edyad(0,0,0,0,0,0,0,0,0,ansFrame); else Eresult := Evector(0,0,0,ansFrame); fi; for Ev in EvecList do Eresult := EaddMixed(1,Eresult,1,Ev,ansFrame) od; RETURN(Eresult) end: ############################################################################## # Just adds components with no accounting for frame ############################################################################## Eadd:= proc(s1,ev1,s2,ev2) [matadd(ev1[1],ev2[1],s1,s2),ev1[2]] end: ############################################################################## # Dyadic operator form of Evector addition ############################################################################## `&++` := proc(ev1,ev2) EaddMixed(1,ev1,1,ev2,ev2[2]) end: ############################################################################## # Dyadic operator form of Evector subtraction ############################################################################## `&--` := proc(ev1,ev2) EaddMixed(1,ev1,-1,ev2,ev2[2]) end: ############################################################################## # multiply an Evector by a scalar ############################################################################## Esm:= proc(s,ev) Eadd(s,ev,0,ev); end: ############################################################################## # dyadic operator form of multiplying an Evector by scalar ############################################################################## `&**` := proc(scal,ev) RETURN(Esm(scal,ev)); end: ############################################################################## # Cross product of two Evectors in same frame ############################################################################## Ecross := proc(ej1,ej2) local ej3; if type([args],[STEvector,STEvector]) then ej3:=express(ej1,ej2[2]); RETURN([crossprod(ej3[1],ej2[1]),ej2[2]]); elif type([args],[STEdyad,STEvector]) then ej3:=express(ej1,ej2[2]); RETURN([multiply(ej3[1],VtoD(ej2)[1]),ej2[2]]); elif type([args],[STEvector,STEdyad]) then ej3:=express(ej1,ej2[2]); RETURN([multiply(VtoD(ej3)[1],ej2[1]),ej2[2]]); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # outer product added as a dirty fix Eoutp := proc(ev1,ev2) local ec1,ec2; if type([args],[STEvector,STEvector]) then ec1:=ev1[1]; ec2:=express(ev2,ev1[2])[1]; RETURN([multiply(convert(ec1,matrix),transpose(convert(ec2,matrix))),ev1[2]]); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # Cross product with Evectors in different frames ############################################################################## EcrossMixed := proc(evec1,evec2,ansFrame) local ev1,ev2; ev1 := express(evec1,ansFrame); ev2 := express(evec2,ansFrame); RETURN(Ecross(ev1,ev2)); end: ############################################################################## # Dyadic operator form of cross product operation ############################################################################## `&xx` := proc(ev1,ev2) #EcrossMixed(ev1,ev2,ev2[2]) RETURN(Ecross(ev1,ev2)); end: ############################################################################## # The magnitude of an Evector ############################################################################## Emag:=proc(Ev) sqrt(Edot(Ev,Ev)); end: ############################################################################## # The antisymmetric Dyad associated with an Evector ############################################################################## VtoD := proc(ev) local vp; vp := &vPart ev; Edyad(0,-vp[3],vp[2],vp[3],0,-vp[1],-vp[2],vp[1],0,&fPart ev); end: ############################################################################## `&VtoD` := proc(ev) VtoD(ev) end: `&DtoV` := proc(ed) DtoV(ed) end: ############################################################################## # The vector associated with the antisym. part of a Dyad ############################################################################## DtoV := proc(ed) local dp,dm; dp := &vPart ed; dp := matadd(dp,transpose(dp),1/2,-1/2); Evector(dp[3,2],dp[1,3],dp[2,1],&fPart ed); end: ############################################################################## # THE CALCULUS OF EVECTORS ############################################################################## ############################################################################## # Sets up the set of substitutions for time differentiation # to have a simple form of result i.e. qt,qtt for diff(q(t),t) ############################################################################## declare:= proc(var) global toTimeFunction, toTimeExpression, varList; if not assigned(toTimeFunction) then toTimeFunction:={} fi; if not assigned(toTimeExpression) then toTimeExpression:={} fi; if not assigned(varList) then varList:={} fi; varList:={op(varList),var}; toTimeFunction := {op(toTimeFunction),var=var(t),cat(var,t)=diff(var(t),t),cat(var,tt)= diff(var(t),t,t)}; toTimeExpression := {op(toTimeExpression),var(t)=var,diff(var(t),t)=cat(var,t),diff(var(t),t,t)=cat(var,tt)}; print(var,`declared`); end: ############################################################################## # ############################################################################## declareList:=proc(varList) local var; for var in varList do declare(var); od; end: ############################################################################## # ############################################################################## dependsTime := proc() local var; for var in [args] do declare(var); od; end: ############################################################################## # ############################################################################## diffTime := proc(expr) local expr1; if (type(expr,algebraic) or type(expr,algebraic=algebraic)) then expr1:=subs(toTimeFunction,expr); expr1:=diff(expr1,t); RETURN(subs(toTimeExpression,expr1)); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # ############################################################################## diffVector := proc(evec,var) local components; if type([args],[STEvector,STuaname]) then components:=SophiaSelectEvectorComponents(evec); RETURN(SophiaConstructEvector(map(diff,components,var),SophiaSelectEvectorFrame(evec))); elif type([args],[STEdyad,STuaname]) then RETURN([map(diff,evec[1],var),evec[2]]); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # ############################################################################## diffVectorTime := proc(evec) local evec1; if type(evec,STEvector) then evec1:=subs(toTimeFunction,evec); evec1:=diffVector(evec1,t); RETURN(subs(toTimeExpression,evec1)); elif type(evec,STEdyad) then evec1:=subs(toTimeFunction,evec); evec1:=diffVector(evec1,t); RETURN(subs(toTimeExpression,evec1)); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # ############################################################################## diffFrame := proc(evec,var,frame) local vFrame,vtran,evectemp; vFrame:=SophiaSelectEvectorFrame(evec); if vFrame=frame then RETURN([map(diff,evec[1],var),frame]) else vtran:=express(evec,frame); evectemp:=[map(diff,vtran[1],var),frame]; RETURN(express(evectemp,vFrame)); fi; end: ############################################################################## # ############################################################################## diffFrameTime := proc(evec,frame) local evec1; if type([args],[STEvector,STFrameName]) then evec1:=express(evec,frame); evec1:=diffVectorTime(evec1); RETURN(express(evec1,SophiaSelectEvectorFrame(evec))); elif type([args],[STEdyad,STFrameName]) then evec1:=express(evec,frame); evec1:=diffVectorTime(evec1); RETURN(express(evec1,evec[2])); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # ############################################################################## `&fdt` := proc(Fr,ev) diffFrameTime(ev,Fr); end: ############################################################################## # ############################################################################## `&dt` := proc(expr) diffTime(expr); end: ############################################################################## # KINEMATIC RELATIONS ############################################################################## genVel := proc(degFrd,qq,uu) local ii,kdee; kdee:={seq(cat(qq,ii,t) = cat(uu,ii), ii=1..degFrd)} end: ############################################################################## # ############################################################################## `&kde` := proc(n) global kde; local ii; dependsTime(seq(cat(q,ii),ii=1..n)); dependsTime(seq(cat(u,ii),ii=1..n)); kde := genVel(n,q,u); print(`Generalized Velocity Form of Kinematic d.e.`); print(` Saved as the GLOBAL Varable kde`); print(kde); end: ############################################################################## # ############################################################################## darBouxVector := proc(A,B,p) local b1,b2,b3,b1d,b2d,b3d; b1 := B&>1; b2 := B&>2; b3 := B&>3; b1d := &simp diffFrame(b1,p,A); b2d := &simp diffFrame(b2,p,A); b3d := &simp diffFrame(b3,p,A); RETURN(Esimplify(B &ev [b2d &o b3, b3d &o b1, b1d &o b2])); end: ############################################################################## # ############################################################################## angularVelocity := proc(B,A) local b1,b2,b3,b1d,b2d,b3d; b1 := Evector(1,0,0,B); b2 := Evector(0,1,0,B); b3 := Evector(0,0,1,B); b1d := diffFrameTime(b1,A); b2d := diffFrameTime(b2,A); b3d := diffFrameTime(b3,A); RETURN(Esimplify(Evector(Edot(b2d,b3),Edot(b3d,b1),Edot(b1d,b2),B))); end: ############################################################################## # ############################################################################## angularAcceleration := proc(B,A,wBA) if wBA[2] = B or wBA[2] = A then RETURN(Esimplify(diffVectorTime(wBA))); else RETURN(Esimplify(diffFrameTime(wBA,A))); fi; end: ############################################################################## # ############################################################################## DiffAngVelOp := proc(evec,AngVel,ansFrame) local evc,evt,evw; evc := Esimplify(express(evec,AngVel[2])); evt := diffVectorTime(evc); evw := EcrossMixed(AngVel,evec,ansFrame); EaddMixed(1,evt,1,evw,ansFrame); end: ############################################################################## # ############################################################################## `&fdtAV` := proc(omega,ev) ((&fPart omega) &fdt ev ) &++ (omega &xx ev) end: ############################################################################## # ############################################################################## aV := proc(fA,fB) local b1,b2,b3,b1d,b2d,b3d; if type([args],[STFrameName,STFrameName]) then b1 := Evector(1,0,0,fB); b2 := Evector(0,1,0,fB); b3 := Evector(0,0,1,fB); b1d := diffFrameTime(b1,fA); b2d := diffFrameTime(b2,fA); b3d := diffFrameTime(b3,fA); RETURN(Esimplify(Evector(Edot(b2d,b3),Edot(b3d,b1),Edot(b1d,b2),fB))); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # ############################################################################## `&aV` := proc(A,B) aV(A,B) end: ############################################################################## # ############################################################################## aA := proc() local fA,fB,omega; fA := args[1]; fB := args[2]; if nargs=3 then omega:=args[3]; else omega:=(fA &aV fB); fi; if (((&fPart omega) = fB) or ((&fPart omega) = fA)) then RETURN(Esimplify(diffVectorTime(omega))); else RETURN(Esimplify(diffFrameTime(omega,fA))); fi; end: ############################################################################## # ############################################################################## FFV := proc(N,rB,wNB,rQ) Esimplify((N &fdt rB) &++ ( wNB &xx (rQ &-- rB))); end: ############################################################################## # ############################################################################## `&ffv` := proc(NB,rBrQ) # call in form [N,B] &NVBQ [rB,rQ] local NN,BB,rBB,rQQ,ww; NN:=NB[1]; BB:=NB[2]; rBB:=rBrQ[1]; rQQ:=rBrQ[2]; ww:= (NN &aV BB); Esimplify((NN &fdt rBB) &++ ( ww &xx (rQQ &-- rBB))); end: ############################################################################## # SCREW ALGEBRA ############################################################################## screwOf:=proc(rvLst,fN) local rivi, ri,vi,vs,screw; vs:=Evector(0,0,0,fN);screw:=Evector(0,0,0,fN); for rivi in rvLst do ri:=rivi[1];vi:=rivi[2]; vs:=Esimplify(EaddList([vi,vs],fN)); screw:=Esimplify(EaddList([EcrossMixed(ri,vi,fN),screw],fN)); od; [screw,vs]; end: ############################################################################## # ############################################################################## screwShift:=proc(screw,eV,rOldNew,fN) Esimplify( EaddMixed(1,screw,1,EcrossMixed(eV,rOldNew,fN),fN)) end: ############################################################################## # ############################################################################## wrench:=proc(screw,eV,fN) local newScrew,rOldNew,er,mag; mag:=simplify(1/Emag(eV)); er:=Esimplify(Esm(mag,eV)); newScrew:=Esimplify(Esm(Edot(screw,er),er)); rOldNew:=Esimplify(Esm(mag,EcrossMixed(er,screw,fN))); [newScrew,rOldNew]; end: ############################################################################## # ############################################################################## `&ss` := proc(screw,shift) local TP,RS; TP := screw[1]; RS := screw[2]; &simp ( TP &++ (RS &xx shift)); end: ############################################################################## # ############################################################################## `&screw` := proc(fN,rvLst) screwOf(rvLst,fN); end: ############################################################################## # DYNAMIC QUANTITIES (SOPHIA 1 STYLE) ############################################################################## ############################################################################## # ############################################################################## linInertiaForce := proc(mass,vel,inertFrame) diffFrameTime(Esm(mass,vel),inertFrame); end: ############################################################################## # ############################################################################## angularMomentum := proc(wBA,IB,B) local w,IBM; w:=B &to wBA; IBM:=B &to IB; Esimplify(IBM &o w); end: ############################################################################## # ############################################################################## rotInertiaForce := proc(wBA,IB,B,A) diffVectorTime(express(angularMomentum(wBA,IB,B),A)); end: ############################################################################## # ############################################################################## rotInertiaForceFrame := proc(wBA,IB,B,A,ansFrame) diffVectorTime(express(angularMomentum(wBA,IB,B),ansFrame)); end: ############################################################################## # BODY PROPERTIES SUCH AS MASS CENTER, MOMENT OF INERTIA ############################################################################## ############################################################################## # ############################################################################## InertiaSet := proc(body,i11,i22,i33,i12,i13,i23) IN||body := array(symmetric,1..3,1..3); IN||body[1,1]:=i11; IN||body[2,2]:=i22; IN||body[3,3]:=i33; IN||body[1,2]:=i12; IN||body[1,3]:=i13; IN||body[2,3]:=i23; print(`declared inertia tensor `, evaln(IN||body)) end: ############################################################################## # ############################################################################## EinertiaDyad := proc(i11,i22,i33,i12,i13,i23,fN) local INbody; if type([args],[algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,STuaname]) then INbody := array(symmetric,1..3,1..3); INbody[1,1]:=i11; INbody[2,2]:=i22; INbody[3,3]:=i33; INbody[1,2]:=i12; INbody[1,3]:=i13; INbody[2,3]:=i23; RETURN([op(INbody),fN]); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # ############################################################################## massPointInertia:=proc(m,rs,fN) local x,y,z,rsN; rsN:=Esimplify(express(rs,fN)); x:=rsN[1][1]; y:=rsN[1][2]; z:=rsN[1][3]; EinertiaDyad(m*(y^2+z^2),m*(x^2+z^2),m*(x^2+y^2),-m*x*y,-m*x*z, -m*y*z,fN); end: ############################################################################## # ############################################################################## massPointSystemInertia:=proc(mPlst,fN) local miri, mi,ri,iner; iner:=EinertiaDyad(0,0,0,0,0,0,fN); for miri in mPlst do mi:=miri[1];ri:=miri[2]; iner:=EaddMixed(1,iner,1,massPointInertia(mi,ri,fN),fN); od; end: ############################################################################## # ############################################################################## InertiaDyadRectangularSolid:=proc(m,xl,yl,zl,fB) local I1,I2,I3; I1:=m*(yl^2+zl^2)/12; I2:=m*(xl^2+zl^2)/12; I3:=m*(xl^2+yl^2)/12; print(`center of mass is at the geometric center`); EinertiaDyad(I1,I2,I3,0,0,0,fB); end: ############################################################################## # ############################################################################## InertiaDyadRightCircularCylinder:=proc(m,R,h,symaxis,fB) local Ir,Is,Ei; Ir:= m*(3*R^2+h^2)/12; Is:= m*R^2/2; print(`center of mass at `, h/2 , ` from bottom or top and on axis`); Ei:=EinertiaDyad(Ir,Ir,Ir,0,0,0,fB); Ei[1][symaxis,symaxis]:=Is; Ei; end: ############################################################################## # ############################################################################## InertiaDyadRightTriangle:=proc(m,base,height,fB) local I1,I2,I12,I3; I1:=m*height^2/18; I2:=m*base^2/18; I12:=m*base*height/36; I3:=I1+I2; print(`center of mass at `,base/3,height/3,` from right corner`); EinertiaDyad(I1,I2,I3,I12,0,0,fB); end: ############################################################################## # ############################################################################## InertiaDyadRectangle:=proc(m,base,baxis,height,haxis,outaxis,fB) local I1,I2,In; I1:=m*height^2/12; I2:=m*base^2/12; print(`center of mass at `,base/2,height/2,` from corner`); In:=EinertiaDyad(I1,I2,I1+I2,0,0,0,fB); In[1][baxis,baxis]:=I1; In[1][haxis,haxis]:=I2; In[1][outaxis,outaxis]:=I1+I2; In; end: ############################################################################## # ############################################################################## `&mpI` := proc(m,rs) massPointInertia(m,rs,&fPart rs) end: ############################################################################## # ############################################################################## `&msI` := proc(frm,mplist) massPointSystemInertia(mplist,frm); end: ############################################################################## # ############################################################################## inertiaAxisTransfer:=proc(istar,ms,rpstar,fN) local imass; imass:=massPointInertia(ms,rpstar,fN); EaddMixed(1,imass,1,istar,fN); end: ############################################################################## # ############################################################################## `&->` := proc(mI,rs) inertiaAxisTransfer(mI[2],mI[1],((-1)&**rs),&fPart rs); end: ############################################################################## # ############################################################################## `&<-` := proc(mIQ,rc) mIQ[2] &-- (mIQ[1] &mpI rc); end: ############################################################################## # ############################################################################## massPointCenter:= proc(massPairList,N) local rs,ms,ex,ri,mi; ms:=0; rs:=Evector(0,0,0,N); for ex in massPairList do ri:= express(ex[2],N); mi:= ex[1]; rs:=Eadd(1,rs,mi,ri); ms:=ms+mi; od; rs:=[Esimplify(Esm(1/ms,rs)),ms]; end: ############################################################################## # ############################################################################## `&cm` := proc(N,mpl) massPointCenter(mpl,N)[1]; end: ############################################################################## # K VECTORS SOPHIA 2 STYLE CALLED KM VECTORS ############################################################################## ############################################################################## # CONSTRUCTORS AND SELECTORS ############################################################################## # Constructs KMvectors ############################################################################## KMvec := proc() RETURN([args,nargs]); end: `&KM` := proc(listOfEvec) KMvec(op(listOfEvec)); end: ############################################################################## # ############################################################################## `&sPartKM` := proc(KMv) RETURN((nops(KMv)-1)); end: ############################################################################## # ############################################################################## `&vPartKM` := proc(KMv) local jjj; [seq(KMv[jjj],jjj=1..(&sPartKM KMv))] end: ############################################################################## # KM ALGEBRA AND SIMPLIFICATION ############################################################################## `&O` := proc(KM1,KM2) local jjj; if (&sPartKM KM1)=(&sPartKM KM2) then RETURN( sum( ( 'KM1[jjj] &o KM2[jjj]'), 'jjj'=1..(&sPartKM KM1) )); else ERROR(`KM size not equal or arg is not a KM vector in &O`); fi; end: ############################################################################## # ############################################################################## KMadd := proc(s1,KM1,s2,KM2) local KM1v,KM2v,KM3v,kk,jjj,ev1,ev2 ; KM1v := (&vPartKM KM1); KM2v := (&vPartKM KM2); KM3v :=[]; kk :=( &sPartKM KM1); for jjj from 1 to kk do ev1 := KM1v[jjj]; ev2 := KM2v[jjj]; KM3v:= [op(KM3v), ((s1 &** ev1) &++ (s2 &** ev2)) ] od; &KM KM3v end: ############################################################################## # ############################################################################## `&+++` := proc(vK1,vK2) if vK1 = [ ] then RETURN(vK2) fi; if VK2 = [ ] then RETURN(vK1) fi; KMadd(1,vK1,1,vK2) end: ############################################################################## # ############################################################################## `&---` := proc(vK1,vK2) KMadd(1,vK1,-1,vK2) end: ############################################################################## # ############################################################################## KMsm := proc(sm,KMvv) local KMv,sz,KMo,jj ; KMo := []; KMv:= &vPartKM KMvv; sz := &sPartKM KMvv; for jj from 1 to sz do KMo :=[op(KMo), (sm &** KMv[jj])] od; &KM KMo end: ############################################################################## # ############################################################################## `&***` := proc(sm,KMvv) KMsm(sm,KMvv) end: ############################################################################## # ############################################################################## `&Ksimp` := proc(KMvv) local KMv,KMo,jj,sz ; KMo := []; KMv:= &vPartKM KMvv; sz := &sPartKM KMvv; for jj from 1 to sz do KMo :=[ op(KMo), (&simp KMv[jj]) ] od; &KM KMo end: KMsimplify := proc(kmv) local n; n:=nops(kmv)-1; KMvec(op(map(Esimplify,&vPartKM kmv))); end: ############################################################################## # ############################################################################## KMtangents:= proc(KMv,var,dimen) local ii,jj,n,imx,zs,KMtvt,vvt,KMtv,eachc,eachev; KMtv:=[]; # argument check not totally good, KMvectors are hard to typecheck if type([args],[list({STEvector,integer}),STuaname,integer]) then zs:={seq(cat(var,ii)=0,ii=1..dimen)}; KMtvt:=subs(zs,KMv); print(`Explicitly time dependent part:`,KMtvt); imx:=array(1..dimen,1..dimen,identity); for ii from 1 to dimen do zs:={seq(cat(var,jj)=imx[ii,jj],jj=1..dimen)}; KMtv:=[op(KMtv),KMadd(1,subs(zs,KMv),-1,KMtvt)]; od; elif type([args],[list({STEvector,integer}),list(STuaname)]) then n:=nops(var); zs:={seq(var[ii]=0,ii=1..n)}; KMtvt:=subs(zs,KMv); print(`Explicitly time dependent part:`,KMtvt); imx:=array(1..n,1..n,identity); for ii from 1 to n do zs:={seq(var[jj]=imx[ii,jj],jj=1..n)}; KMtv:=[op(KMtv),KMadd(1,subs(zs,KMv),-1,KMtvt)]; od; else ERROR(`Invalid arguments!`); fi; RETURN(KMtv); end: ############################################################################## # ############################################################################## vKtime := proc(KMv,var,dimen) local zs,n,ii; # argument check not totally good, KMvectors are hard to typecheck if type([args],[list({STEvector,integer}),STuaname,integer]) then zs:={seq(cat(var,ii)=0,ii=1..dimen)}; RETURN(subs(zs,KMv)); elif type([args],[list({STEvector,integer}),list(STuaname)]) then n:=nops(var); zs:={seq(var[ii]=0,ii=1..n)}; RETURN(subs(zs,KMv)); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # KM VECTOR CALCULUS ############################################################################## `&Kfdt` := proc(frm,KMvv) local KMv,KMo,jj,sz ; KMo := []; KMv:= &vPartKM KMvv; sz := &sPartKM KMvv; for jj from 1 to sz do KMo :=[ op(KMo), (frm &fdt KMv[jj]) ] od; &KM KMo; end: ############################################################################## # ############################################################################## `&rKfdt` := proc(frm,KMvv) local KMv,KMvE,KMo,jj,sz ; KMo := []; KMv:= &vPartKM KMvv; sz := &sPartKM KMvv; for jj from 1 to sz do KMvE := KMv[jj]; if type(KMvE,symbol) then KMo := [op(KMo),frm &aV KMvE]; else KMo :=[ op(KMo), (frm &fdt KMv[jj]) ]; fi; od; &KM KMo end: ############################################################################## # KM Mass Product # makes momentum from mass and velocity ############################################################################## KMMP:=proc(m,v) local n,ii,w; n:=nops(v)-1; w:=[]; for ii from 1 to n do w:=[op(w),Edot(m[ii],v[ii])]; od; w:=[op(w),n]; end: ############################################################################## # KMip # inner product on the tangent vectors in the mass metric ############################################################################## KMIP:=proc(u,v,m) local n,ii,w; n:=nops(u)-1; if (n=v[nops(v)] and n=m[nops(m)]) then w:=KMMP(m,v); RETURN(sum(('u[ii] &o w[ii]'),ii=1..n)); else ERROR(`Invalid arguments!`) fi; end: ############################################################################## # KM Scalar Product ############################################################################## KMSP:=proc(s,v) local n,ii,w; n:=nops(v)-1; w:=[]; for ii from 1 to n do w:=[op(w),Esm(s,v[ii])]; od; w:=[op(w),n]; end: ############################################################################## # Mass Tangent Product # Takes mass operator and list of tangent KMvectors ############################################################################## MBP:=proc(m,b) local n,ii,p; n:=nops(b); p:=[]; for ii from 1 to n do p:=[op(p),KMMP(m,b[ii])]; od; RETURN(p); end: ############################################################################## # SUPER K VECTORS (LISTS OF K VECTORS) ############################################################################## `&SK` := proc(KMlist) KMlist end: ############################################################################## # ############################################################################## `&SKc` := proc(Vsk,j) Vsk[j] end: ############################################################################## # ############################################################################## `&SKn` := proc(Vsk) nops(Vsk) end: ############################################################################## # ############################################################################## `&****` := proc(sc,Vsk) local Wsk,Vk,jj,nn; Wsk:= [NULL]; nn := ( &SKn Vsk); for jj from 1 to nn do Vk := (Vsk &SKc jj); Wsk:= [op(Wsk),(sc &*** Vk)] od; &SK Wsk end: ############################################################################## # ############################################################################## `&++++` := proc(Vsk1,Vsk2) local Wsk, Vk1, Vk2, nn,jj; if Vsk1=[ ] then RETURN(Vsk2) fi; if Vsk2=[ ] then RETURN(Vsk1) fi; Wsk := [NULL]; nn := (&SKn Vsk1); for jj from 1 to nn do Vk1 := (Vsk1 &SKc jj); Vk2 := (Vsk2 &SKc jj); Wsk := [op(Wsk),(Vk1 &+++ Vk2)] od; &SK Wsk end: ############################################################################## # ############################################################################## `&----` := proc(Vsk1,Vsk2) Vsk1 &++++ ( (-1) &**** Vsk2 ) end: ############################################################################## # ############################################################################## `&++**` := proc(BetaTau,Tau) local jj,ii,nn,Beta,betaii; nn := &SKn Tau; Beta := [ ]; for ii from 1 to nn do betaii := [ ]; for jj from 1 to nn do betaii:= ( betaii &+++ (BetaTau[ii,jj] &*** Tau[jj]) ); od; Beta := [op(Beta),betaii]; od; Beta; end: ############################################################################## # METRICS AND TRANSFORMATIONS OF K VECTORS ############################################################################## `&Kmetric`:=proc(tauK) local Gtau,nz,ii,jj; nz := &SKn tauK; Gtau := array(symmetric,1..nz,1..nz); for ii from 1 to nz do for jj from ii to nz do Gtau[ii,jj]:= tauK[ii] &O tauK[jj]; od; od; op(Gtau); end: ############################################################################## # ############################################################################## `&cv` := proc(tauK) local Gt,GtC,tC; Gt := &Kmetric tauK; GtC := inverse(Gt); tC := op(GtC) &++** tauK; end: ############################################################################## # ############################################################################## `&<>4` := proc(ctauK,wwK) # coef of wwK in set tauK # note that difference from &kane is that cf is an array local nz,cf,ii; nz := &SKn ctauK; cf := array(1..nz); for ii from 1 to nz do cf[ii]:= ctauK[ii] &O wwK od; op(cf); end: ############################################################################## # ############################################################################## `&SKsum` := proc(cf,tauK) local nz,sK,ii; nz := &SKn tauK; sK := cf[1] &*** tauK[1]; for ii from 2 to nz do sK := sK &+++ (cf[ii] &*** tauK[ii]) od; sK end: ############################################################################## # ############################################################################## `&kane` := proc(tau,wK) local tauj,tw,eqst; eqst := [ ]; for tauj in tau do tw := tauj &O wK; eqst := [op(eqst),tw]; od; RETURN(eqst); end: ############################################################################## # UTILITY FUNCTIONS ############################################################################## clearVars := proc() global toTimeFunction, toTimeExpression, varExpression, varFunction, varList; toTimeFunction := 'toTimeFunction'; toTimeExpression := 'toTimeExpression'; varExpression := 'varExpression'; varFunction := 'varFunction'; varList := 'varList'; lprint(`varFunction, varExpression, toTimeFunction, toTimeExpression cleared`); lprint(``); end: ############################################################################## # ############################################################################## plist:=proc(state,j,k,t1,t2,dt) local ti,st,ls; ls:= NULL; for ti from eval(t1) by eval(dt) to eval(t2) do st:=state(ti); ls:=ls,st[j],st(ti)[k] od; ls:=[ls]; end: ############################################################################## # ############################################################################## setSymbolicTrig:=proc(var,numb) global toTrigSymbolic, toTrigExplicit; local i; toTrigSymbolic:={}; toTrigExplicit:={}; for i from 1 to numb do toTrigSymbolic:= toTrigSymbolic union {sin(cat(var,i))=cat(s,i),cos(cat(var,i))=cat(c,i)}; toTrigExplicit:= toTrigExplicit union {cat(s,i)=sin(cat(var,i)),cat(c,i)=cos(cat(var,i))}; od end: ############################################################## # Maple-Sophia Extension for Cartesian Based Angular Velocity Computations # M.B. Lesser, Nov 1993 # KTH-Stockholm ########################################################################## # rb1, rb2 and rb3 are given Evector positions, generally functions of time # also being assummed that they go to points fixed in the body and frA is # the frame in which the reference point is given, frB the body frame. # with 5 arguments it places the defined rotation matrix and frames on # the list of rotation sets, with 3 arguments it simply outputs the # triad fixed in the body frame. # # Hacked for SophiaConstruct... etc ############################################################################ PointSetTriad := proc(rb1,rb2,rb3,frA,frB) global SophiaGlobalDeclaredFrameRelations; local zb1,zb2,zb3,b1,b2,b3; zb1 := rb2 &-- rb1; b1 := (1/Emag(zb1)) &** zb1; zb2 := rb3 &-- rb1; zb3 := zb2 &-- ((zb2 &o b1) &** b1); b2 := (1/Emag(zb3)) &** zb3; b3 := b1 &xx b2; if nargs = 5 then b1 := frA &to b1; b2 := frA &to b2; b3 := frA &to b3; SophiaConstructAddRmx(frA,frB,b1 &c 1,b1 &c 2,b1 &c 3,b2 &c 1,b2 &c 2,b2 &c 3,b3 &c 1,b3 &c 2,b3 &c 3); SophiaGlobalDeclaredFrameRelations:=SophiaGlobalDeclaredFrameRelations union {{frA,frB}}; fi; RETURN([b1,b2,b3]); end: `&PST` := proc(frs,pts) PointSetTriad(pts[1],pts[2],pts[3],frs[1],frs[2]); end: ############################################################################## # Sets up rotations in terms of all the direction cosine # matrix elements. The parameter name prefix is assigned by # the user as the third argument. # the forth argument is the angular velocity vector wAB as a # list of components (not an Evector) # as represented in either frame A or B ############################################################################## poissonRot := proc(A,B,s,av,avFrame) global SophiaGlobalDeclaredFrameRelations; local omega,poisson,poissonSet,i,j; SophiaConstructAddRmx(B,A,cat(s,11),cat(s,12),cat(s,13),cat(s,21),cat(s,22),cat(s,23),cat(s,31),cat(s,32),cat(s,33)); SophiaGlobalDeclaredFrameRelations:=SophiaGlobalDeclaredFrameRelations union {{A,B}}; declareList([cat(s,11),cat(s,12),cat(s,13),cat(s,21),cat(s,22),cat(s,23),cat(s,31),cat(s,32),cat(s,33)]); omega:= matrix([[0,-av[3],av[2]],[av[3],0,-av[1]],[-av[2],av[1],0]]); if avFrame=A then poisson:=map(diffTime,Rmx(A,B))=evalm(omega &* Rmx(A,B)) fi; if avFrame=B then poisson:=map(diffTime,Rmx(A,B))=evalm(Rmx(A,B) &* omega) fi; poissonSet := { }; for i from 1 to 3 do for j from 1 to 3 do poissonSet:=poissonSet union {lhs(poisson)[i,j]=rhs(poisson)[i,j]}; od; od; RETURN(poissonSet); end: ############################################################################## # ############################################################################## dependsVar:= proc(var,p) global varFunction, varExpression, varList; if not assigned(varFunction) then varFunction:={} fi; if not assigned(varExpression) then varExpression:={} fi; if not assigned(varList) then varList:={} fi; varList:={op(varList),var}; varFunction :={op(varFunction),var=var(p),cat(var,p)=diff(var(p),p),cat(var,p,p)= diff(var(p),p,p)}; varExpression :={op(varExpression),var(p)=var,diff(var(p),p)=cat(var,p),diff(var(p),p,p)=cat(var,p,p)}; lprint(var,`declared`) end: ############################################################################## # ############################################################################## `&dependsOn` := proc(var,p) dependsVar(var,p) end: ############################################################################## # ############################################################################## diffVar := proc(expr,var) local expr1; expr1:=subs(varFunction,expr); expr1:=diff(expr1,var); subs(varExpression,expr1) end: ############################################################################## # ############################################################################## `&dVar` := proc(expr,var) diffVar(expr,var) end: ############################################################################## # ############################################################################## CoefficientArray:=proc(exprList,termList) local expr1,cfArray,ii,jj,nr,nc,trm; nr:= nops(exprList); nc:= nops(termList); cfArray := matrix(nr,nc); for ii from 1 to nr do for jj from 1 to nc do trm := termList[jj]; cfArray[ii,jj]:= coeff(collect(exprList[ii],termList),trm); od; od; eval(cfArray); end: ############################################################################## # ############################################################################## SourceList:= proc(exprList,termList) local zList,trm; zList:=[ ]; for trm in termList do zList:= [op(zList),trm=0] od; subs(zList,exprList); end: ############################################################################## # Sophia Routine to Produce the Reduced Kinematic Differential Equations # For a system that is either Nonholonomic or Redundant ############################################################################## # © Martin Lesser October 12, 1994 ############################################################################## # Usage: ############################################################################## # arguments # kde is a set of kinematic differential equations # vce is a set of velocity constraint relations # gs is a LIST of the basic group of generalized speeds that are to be # retained. This must equal the independent degrees of freedom # us is the name of the generalized speed symbol, usually u # ws is the name for the transformed generalized speeds ############################################################################## # example # pendulum of length L using cartesian coordinates q1, q2 # kde := {q1t=u1,q2t=u2} # vce := { q1t * q1- (L-q2)*q2t = 0 which can also be given in terms of u1,u2 # gs := [u1] # us := u # ws := w ############################################################################## # The result will be a set of kinematic differential equations # for q1t, q2t in terms of a single generalized speed w1, thus # {q2t = q1*w1/(L-q2), q1t = w1} ############################################################################## # The general use of this is to use the velocity constraint relations to # obtain the correct number of tangent Kvectors for the problem. # It is also used as a substitution in all differentiations to eliminate qt terms. # If ng is the geometric degree of freedom, i.e. the number of kde # nc the number of velocity constraint equations # ni = ng - nc # and there must be ni generalized speeds in bgs # THIS is checked! if not so an error message is given # The user must be careful that bgs contains variables which are # independent, if not nonsense can result which may or may not # show up as an error in the rest of the procedure. ############################################################################## # Modifications History # Oct 18, 1994 version 0.91 # new function ReduceGSKDE output to a sequence of sets # rigst, rtkde # that is the reduced inverse generalized speed transformation and # the reduced transformed kinematic differential equations # changed nameing convention bgs to bg and dce to vce ############################################################################## kdeReduce := proc(kde,vce,gs,us,ws) local ng,nc,ni,ni1,gst,igst,rigst,j; ng := nops(kde); nc := nops(vce); ni := nops(gs); ni1 := ni+1; if not(ng = nc + ni) then ERROR(`inconsistant numbers in kde, vce and bgs`) fi; gst:={seq(cat(ws,j)=gs[j],j=1..ni)} union {seq(cat(ws,j)=subs(kde,lhs(vce[j-ni])),j=ni1..ng )}; igst:= solve(gst,{seq(cat(us,j),j=1..ng)}); rigst:= subs({seq(cat(ws,j)=0,j=ni1..ng)},igst); subs(rigst,kde); end: ReducedSpeeds := proc(gs,kde,vce,us,ws) local ng,nc,ni,ni1,gst,igst,rigst,j; ng := nops(kde); nc := nops(vce); ni := nops(gs); ni1 := ni+1; if not(ng = nc + ni) then ERROR(`inconsistant numbers in kde, vce and gs`) fi; gst:={seq(cat(ws,j)=gs[j],j=1..ni)} union {seq(cat(ws,j)=subs(kde,lhs(vce[j-ni])),j=ni1..ng )}; igst:= solve(gst,{seq(cat(us,j),j=1..ng)}); rigst:= subs({seq(cat(ws,j)=0,j=ni1..ng)},igst); rigst,subs(rigst,kde); end: ############################################################################## # var1eqs # computes the derivatives of the vector field with respect to x, # i.e. the first order variational equations. # the resulting equations are returned in a list, following the # convention in MATLAB with the elements of first column followed # by the second column and so on ############################################################################## var1eqs:=proc(ft,vars) local n,ii,jj,eqm,vmx,mxeq; if type([args],[list(name=algebraic),list(name)]) and nops(ft)=nops(vars) then n:=nops(ft); eqm:=jacobian([seq(rhs(ft[ii]),ii=1..n)],vars); vmx:=matrix(n,n,[seq(seq(cat(`d`,vars[jj],`d`,vars[ii]),ii=1..n),jj=1..n)]); mxeq:=multiply(eqm,vmx); RETURN([seq(seq(cat(`d`,vars[ii],`d`,vars[jj],`_t`)=mxeq[ii,jj],ii=1..n),jj=1..n)]); else ERROR(`Wrong arguments!`); fi; end: ############################################################################## # jaceqs # computes the jacobian of xt with respect to x and returns a # list of equations with these new equations ############################################################################## jaceqs:=proc(ft,vars) local ii,jj,n; if type([args],[list(name=algebraic),list(name)]) and nops(ft)=nops(vars) then n:=nops(ft); RETURN([seq(seq(cat(`d`,lhs(ft[ii]),`d`,vars[jj])=diff(rhs(ft[ii]),vars[jj]),ii=1..n),jj=1..n)]); else ERROR(`Invalid arguments!`); fi; end: ############################################################################## # Sophia Routines to handle the development of approximate equations of # motion using Kane's method. ############################################################################## # © Martin Lesser October 31, 1994 # converted for Sophia II by Anders Lennartsson 1996-11-25 ############################################################################## # The purpose of these routines is to supply Sophia Tools for dealing with # solutions near some known solution, most often this is a case of no motion or # displacement. It also includes complete linearization such as discussed in # Kane and Levinson. ############################################################################## # The first series of procedures simplifies the expansion of Sophia Evectors, # Kvector and Svector objects in terms of a parameter. The calling function in all # cases requries three arguments, the object to be expanded , the expansion parameter # and the order NEGLECTED. It is assumed that a substitution set has been created # to change to the perturbed form. For example: # disturbance := { q1= Q1 + eps* q1, u1 = U1 + eps*u1} # and that Q1, U1 are an exact solution so q1 and u1 become perturbations # Note one could use e.g. x1 and y1 for q1 and u1 on the right side but frequently it is # convenient to retain the same name for the perturbation as the original # quantity. # thus EvectSeries(subs(disturbance,v1),eps,3) will produce the velocity Evector up # to and including terms of order 2. # Truncation of products and so forth must be carried out as needed, the package is not # set up to deal with details, just to provide help in manipulation! ############################################################################## EvectSeries := proc(evect,eps1,ns) (&fPart evect) &ev [ convert(series(evect &c 1,eps1=0,ns),polynom),convert(series(evect &c 2,eps1=0,ns),polynom),convert(series(evect &c 3,eps1=0,ns),polynom)]; end: ############################################################################## KvectSeries := proc(kvect,eps1,ns) local kvect1,kvect2,kk,jjj,ev1; kvect1 := (&vPartKM kvect); kvect2 :=[]; kk :=( &sPartKM kvect); for jjj from 1 to kk do ev1 := kvect1[jjj]; kvect2:= [op(kvect2),EvectSeries(ev1,eps1,ns) ] od; &KM kvect2 end: ############################################################################## SvectSeries := proc(svect,eps1,ns) local svect1,svect2,ja,jb,kv1; svect1 := svect; svect2 := []; ja := &SKn svect; for jb from 1 to ja do kv1:= svect1 &SKc jb; svect2 := [op(svect2),KvectSeries(kv1,eps1,ns)] od; &SK svect2 end: ############################################################################## # Transforming coordinates and objects means that we must simplify the # the rotation matrices. This is also required for frame based differentiation # and you should be sure to declare the time dependance of relavant # new variables! ############################################################################## # The first routine sets the stage for the approximate rotational # changes. The procedure SetApproximateRotation takes two # arguments, the disturbance set above and the order of expansion # from which it sets up two groups of global rotation matrices # for all known transformations. For example RABext is the # full or exact matrix for the RAB transformation and # RABapx is the approximate one # The SWitchExactApproximate puts one or the other into # the symbols RAB etc, where it is used by coordinate # transformations and frame based differentiations. # SwitchExactApproximate takes one argument the # character E typed `E` or A typed `A` for exact # or approximate. # The user need not use it in most cases as it is # incorporated in frame based differentiation operators # for approximate analysis, N &Afdt Evector and N &AKfdt Kvector. # These restore the exact expressions after doing their work. ############################################################################## SetApproximateRotation:=proc(disturbance,eps,ns) global SophiaGlobalRmxTable,SophiaGlobalRmxTableApprox,SophiaGlobalRmxTableExact; local each,incs; SophiaGlobalRmxTableExact:=copy(SophiaGlobalRmxTable); SophiaGlobalRmxTableApprox:=subs(disturbance,copy(SophiaGlobalRmxTable)); incs:=indices(SophiaGlobalRmxTable); for each in incs do SophiaGlobalRmxTableApprox[op(each)]:=map(convert,map(series,SophiaGlobalRmxTableApprox[op(each)],eps,ns),polynom); od; RETURN(); end: ############################################################################## SwitchExactApproximate := proc(kind) global SophiaGlobalRmxTable,SophiaGlobalRmxTableApprox,SophiaGlobalRmxTableExact; if not( (kind = `E`) or (kind = `A`) ) then ERROR(`Invalid arguments, only input string E or A allowed in SwitchExactApproximate!`) elif not(assigned(SophiaGlobalRmxTableApprox)) then ERROR(`Must use SetApproximateRotation before Switch`) elif not(assigned(SophiaGlobalRmxTableExact)) then ERROR(`Must use SetApproximateRotation before Switch`) else if kind = `E` then SophiaGlobalRmxTable:=copy(SophiaGlobalRmxTableExact); elif kind =`A` then SophiaGlobalRmxTable:=copy(SophiaGlobalRmxTableApprox); fi; RETURN(); fi; end: # These are the approximate differentiation routines ############################################################################## `&Afdt` := proc(frn,evec) local devec: SwitchExactApproximate(`A`): devec := N &fdt evec: SwitchExactApproximate(`E`): devec end: ############################################################################## `&AKfdt` := proc(frn,kvec) local dkvec: SwitchExactApproximate(`A`): dkvec := N &Kfdt kvec: SwitchExactApproximate(`E`): dkvec end: ############################################################################## # approximate products, an alternative is to switch overall environment to # approximate status `&Akane` := proc(svect,kvect) local lst: SwitchExactApproximate(`A`): lst := svect &kane kvect: SwitchExactApproximate(`E`): lst end: ############################################################################## # These are the approximate frame transfer expressions `&Ato` := proc(frm,objkt) local obk: SwitchExactApproximate(`A`): obk := frm &to objkt: SwitchExactApproximate(`E`): obk; end: ############################################################################## # global variables ############################################################################## SophiaGlobalDeclaredFrameRelations:={}: SophiaGlobalRmxTable:=table(): SophiaGlobalRmxTableExact:='SophiaGlobalRmxTableExact': SophiaGlobalRmxTableApprox:='SophiaGlobalRmxTableApprox': ############################################################################## # INTERNAL PART # TREE OF DECLARED FRAME RELATIONS USED TO CALCULATE NEW Rmx # # ############################################################################## # SophiaSelectRmx(fN1,fN2) # Returns rotation matrix if it exists or may be calculated. ############################################################################## SophiaSelectRmx := proc(fN1,fN2) global SophiaGlobalRmxTable; local path,nn,ii; if (not(type([args],[STFrameName,STFrameName]))) then ERROR(`Invalid arguments!`); fi; if (SophiaCheckIfRmx(fN1,fN2)) then RETURN(SophiaGlobalRmxTable[fN1,fN2]); else path:=SophiaSelectFindDeclaredFramePath(fN1,fN2); if (nops(path)=0) then ERROR(cat(`No declared path between `,fN1,`and`,fN2,`.`)); fi; nn:=nops(path); for ii from 3 to nn do if not(SophiaCheckIfRmx(fN1,path[ii])) then SophiaGlobalRmxTable[fN1,path[ii]]:=evalm(SophiaGlobalRmxTable[fN1,path[ii-1]] &* SophiaGlobalRmxTable[path[ii-1],path[ii]]); SophiaGlobalRmxTable[path[ii],fN1]:=copy(transpose(eval(SophiaGlobalRmxTable[fN1,path[ii]]))); fi; od; RETURN(SophiaGlobalRmxTable[fN1,fN2]); fi; end: ############################################################################## # SophiaCheckIfRmx(fN1,fN2) # Returns true if a rotation matrix relating the two frames # have been calculated, otherwise false. ############################################################################## SophiaCheckIfRmx:=proc(fN1,fN2) if (not(type([args],[STFrameName,STFrameName]))) then ERROR(`Invalid arguments!`); else RETURN(assigned(SophiaGlobalRmxTable[fN1,fN2])); fi; end: ############################################################################## # SophiaSelectFindDeclaredFramePath(fN1,fN2,f_set) # Returns a list of declared frames to traverse when going from # one frame to another. # The last argument is optional, used during recursive calls. ############################################################################## `type/SophiaSelectFindDeclaredFramePathargs`:={[STFrameName,STFrameName],[STFrameName,STFrameName,STDeclaredFrameRelationSet]}: SophiaSelectFindDeclaredFramePath:=proc(fN1,fN2,f_set) local frame_set,new_frame_set,each,partners,path_list,scrap; if (nargs=2) then if assigned(SophiaGlobalDeclaredFrameRelations) then frame_set:=SophiaGlobalDeclaredFrameRelations; else RETURN([]); fi; elif (nargs=3) then frame_set:=f_set; else ERROR(`Wrong number of arguments!`); fi; if not(SophiaCheckIfFrame(fN1,frame_set) and SophiaCheckIfFrame(fN2,frame_set)) then RETURN([]); elif SophiaCheckIfDeclaredFrameRelation(fN1,fN2,frame_set) then RETURN([fN1,fN2]); fi; partners:=SophiaSelectFindFramePartners(fN1,frame_set); new_frame_set:={}; for each in frame_set do if not(member(fN1,each)) then new_frame_set:=new_frame_set union {each}; fi; od; path_list:=[]; for each in partners do scrap:=SophiaSelectFindDeclaredFramePath(each,fN2,new_frame_set); if not(nops(scrap)=0) then path_list:=[op(path_list),scrap]; fi; od; if nops(path_list)=0 then RETURN([]); elif nops(path_list)=1 then RETURN([fN1,op(path_list[1])]); else ERROR(`Loops detected in frame definitions, reachable via two ways!`); fi; end: ############################################################################## SophiaCheckIfFrame:=proc(fN1,f_set) local each,test1,frame_set; if (nargs=1) then if assigned(SophiaGlobalDeclaredFrameRelations) then frame_set:=SophiaGlobalDeclaredFrameRelations; else RETURN(false); fi; elif (nargs=2) then frame_set:=f_set; else ERROR(`Wrong number of arguments!`); fi; test1:=false; for each in frame_set do if member(fN1,each) then test1:=true; fi; od; RETURN(test1); end: ############################################################################## SophiaCheckIfDeclaredFrameRelation:=proc(fN1,fN2,f_set) local each,test1,frame_set; if (nargs=2) then if assigned(SophiaGlobalDeclaredFrameRelations) then frame_set:=SophiaGlobalDeclaredFrameRelations; else RETURN(false); fi; elif (nargs=3) then frame_set:=f_set; else ERROR(`Wrong number of arguments!`); fi; test1:=false; for each in frame_set do if (member(fN1,each) and member(fN2,each)) then test1:=true; fi; od; RETURN(test1); end: ############################################################################## SophiaSelectFindFramePartners:=proc(fN,f_set) local frame_set,each,list1; if (nargs=1) then if assigned(SophiaGlobalDeclaredFrameRelations) then frame_set:=SophiaGlobalDeclaredFrameRelations; else RETURN([]); fi; elif (nargs=2) then frame_set:=f_set; else ERROR(`Wrong number of arguments!`); fi; list1:=[]; for each in frame_set do if member(fN,each) then list1:=[op(list1),op(each minus {fN})]; fi; od; list1; end: ############################################################################## erasordef:=proc(f1,f2) SophiaConstructEraseRotDef(f1,f2); end: SophiaConstructEraseRotDef:=proc(f1,f2) local ns1,ns2; global SophiaGlobalRmxTable, SophiaGlobalDeclaredFrameRelations; if not(type([args],[STFrameName,STFrameName])) then ERROR(`Invalid arguments!`); fi; ns1:=SophiaSelectFindFramePartners(f1); ns1:=convert(ns1,set) minus {f2}; if nops(ns1)<>1 then ERROR(`First frame argment adjacent to other frames!`); fi; ns2:=SophiaSelectFindFramePartners(f2); ns2:=convert(ns2,set) minus {f1}; if nops(ns2)<>1 then ERROR(`Second frame argment adjacent to other frames!`); fi; if assigned(SophiaGlobalRmxTable[f1,op(ns2)]) then ERROR(`Definition already used in other matrices, restart necessary if definition in error!`); fi; if assigned(SophiaGlobalRmxTable[f2,op(ns1)]) then ERROR(`Definition already used in other matrices, restart necessary if definition in error!`); fi; SophiaGlobalRmxTable[f1,f2]:='SophiaGlobalRmxTable[f1,f2]'; SophiaGlobalRmxTable[f2,f1]:='SophiaGlobalRmxTable[f2,f1]'; SophiaGlobalDeclaredFrameRelations:=SophiaGlobalDeclaredFrameRelations minus {{f1,f2}}; RETURN(); end: ############################################################################## SophiaConstructAddRmx := proc(sf,nf,b1a1,b1a2,b1a3,b2a1,b2a2,b2a3,b3a1,b3a2,b3a3) global SophiaGlobalRmxTable; if evalb(nargs = 11) then SophiaGlobalRmxTable[nf,sf]:=array(1..3,1..3,[[b1a1,b1a2,b1a3],[b2a1,b2a2,b2a3],[b3a1,b3a2,b3a3]]); SophiaGlobalRmxTable[sf,nf]:=array(1..3,1..3,[[b1a1,b2a1,b3a1],[b1a2,b2a2,b3a2],[b1a3,b2a3,b3a3]]); fi; end: ############################################################################## # SophiaConstructCalcRmx(fN,orspc) # Calculates a rotation matrix from the arguments and adds it to the # global table of rotation matrices. # The indices are determined by given frame names. ############################################################################## SophiaConstructCalcRmx := proc(nfN,orspc) local e1,e2,e3,frame1,axis,angle,r11,r12,r13,r21,r22,r23,r31,r32,r33,vec; if (not(type([args],[STFrameName,STOrientationSpecifier]))) then ERROR(`Invalid arguments!`); fi; frame1:=SophiaSelectFrame(orspc); # axis determined by integer if type(orspc,STOrientationSpecifierSimple) then axis:=SophiaSelectAxis(orspc); angle:=SophiaSelectAngle(orspc); if evalb(axis=1) then SophiaConstructAddRmx(frame1,nfN,1,0,0,0,cos(angle),sin(angle),0,-sin(angle),cos(angle)); RETURN(true); elif evalb(axis=2) then SophiaConstructAddRmx(frame1,nfN,cos(angle),0,-sin(angle),0,1,0,sin(angle),0,cos(angle)); RETURN(true); elif evalb(axis=3) then SophiaConstructAddRmx(frame1,nfN,cos(angle),sin(angle),0,-sin(angle),cos(angle),0,0,0,1); RETURN(true); else ERROR(`Axis integer different than 1,2,3.`); fi; # rotation axis determined by components in starting frame elif type(orspc,STOrientationSpecifierAxis) then vec:=SophiaSelectAxis(orspc); vec:=scalarmul(vec,1/sqrt(dotprod(vec,vec,orthogonal))); angle:=SophiaSelectAngle(orspc); e1:=vec[1]; e2:=vec[2]; e3:=vec[3]; r11:=simplify(e1*e1*(1-cos(angle))+cos(angle)); r12:=simplify(e1*e2*(1-cos(angle))+e3*sin(angle)); r13:=simplify(e1*e3*(1-cos(angle))-e2*sin(angle)); r21:=simplify(e2*e1*(1-cos(angle))-e3*sin(angle)); r22:=simplify(e2*e2*(1-cos(angle))+cos(angle)); r23:=simplify(e2*e3*(1-cos(angle))+e1*sin(angle)); r31:=simplify(e3*e1*(1-cos(angle))+e2*sin(angle)); r32:=simplify(e3*e2*(1-cos(angle))-e1*sin(angle)); r33:=simplify(e3*e3*(1-cos(angle))+cos(angle)); SophiaConstructAddRmx(frame1,nfN,r11,r12,r13,r21,r22,r23,r31,r32,r33); RETURN(true); else ERROR(`This type of orientation specifier not implemented!`); fi; end: ############################################################################## # SophiaSelectFrame(orspc) # Selects the frame from a orientation specifier. ############################################################################## SophiaSelectFrame := proc(orspc) if (not(type([args],[STOrientationSpecifier]))) then ERROR(`Invalid arguments!`); else RETURN(orspc[1]); fi; end: ############################################################################## # SophiaSelectAxis(orspc) # Selects the axis from a orientation specifier of type # STOrientationSpecifierAxis # STOrientationSpecifierSimple ############################################################################## SophiaSelectAxis := proc(orspc) if (not(type([args],{[STOrientationSpecifierAxis],[STOrientationSpecifierSimple]}))) then ERROR(`Invalid arguments!`); else RETURN(orspc[2]); fi; end: ############################################################################## # SophiaSelectAngle(orspc) # Selects the anlge or rotation from a orientation specifier of type # STOrientationSpecifierAxis # STOrientationSpecifierSimple ############################################################################## SophiaSelectAngle := proc(orspc) if (not(type([args],{[STOrientationSpecifierAxis],[STOrientationSpecifierSimple]}))) then ERROR(`Invalid arguments!`); else RETURN(orspc[3]); fi; end: ############################################################################## # SophiaDeclareFrameRelation := proc(ospc) # Declares a relation between two frames if # a) none of them are declared earlier # b) only one of them is declared earlier # c) both declared earlier but no path of declared # frames between them exist ############################################################################## SophiaDeclareFrameRelation := proc(nfN,orspc) global SophiaGlobalDeclaredFrameRelations; local sfr,ffr,frame1; if (not(type([args],[STFrameName,STOrientationSpecifier]))) then ERROR(`Invalid arguments!`); fi; frame1:=SophiaSelectFrame(orspc); sfr:=SophiaCheckIfFrame(frame1); ffr:=SophiaCheckIfFrame(nfN); if ((sfr and ffr) and (nops(SophiaSelectFindDeclaredFramePath(frame1,nfN))>0)) then ERROR(`A declared path between frames`,frame1,` and`,nfN,` already exists!`); fi; if (SophiaConstructCalcRmx(nfN,orspc)) then print(cat(`Frame relation between `,frame1,` and `,nfN,` defined!`)); SophiaGlobalDeclaredFrameRelations:=SophiaGlobalDeclaredFrameRelations union {{frame1,nfN}}; RETURN(true); else print(cat(`Warning! Frame relation between `,frame1,` and `,nfN,` NOT defined!`)); RETURN(false); fi; end: ############################################################################## # # TTTTTTT Y Y PPPPPPP EEEEEEE SSSSSSS # T Y Y P P E S # T Y Y P P E S # T Y PPPPPPP EEEEE SSSSSSS # T Y P E S # T Y P E S # T Y P EEEEEEE SSSSSSS # ############################################################################## # User Types ############################################################################## `type/ST3Index`:={identical(1),identical(2),identical(3)}: `type/STComponent`:=algebraic: `type/STParameter`:={numeric,symbol}: `type/STAngle`:=algebraic: `type/STuaname`:=proc(v) evalb( type(v,symbol) and not(evalb(assigned(v)))) end: `type/STFrameName`:=proc(v) evalb( type(v,symbol) and not(evalb(assigned(v)))) end: `type/STVector`:=proc(v,n) evalb( type(v,'vector'(algebraic)) and evalb(vectdim(v)=n)) end: `type/STSquareMatrix` := proc(mx,n) evalb( type(mx,'matrix'(algebraic,square)) and evalb(rowdim(mx)=n)) end: `type/ST3Vector`:=STVector(3): `type/ST3Dyad`:=STSquareMatrix(3): `type/ST3VList`:=[algebraic,algebraic,algebraic]: `type/ST3DList`:=[algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,algebraic,algebraic]: ############################################################################## # Sophia Types ############################################################################## `type/STEvector`:=[ST3Vector,STFrameName]: `type/STEdyad`:=[ST3Dyad,STFrameName]: `type/STKvector`:=list(STEvector): `type/STAxisInteger`:={identical(1),identical(2),identical(3)}: `type/STAxis`:={STAxisInteger,STAxisEvector,ST,STAxisComponents}: `type/STDeclaredFrameRelation`:=set(STFrameName,STFrameName): `type/STDeclaredFrameRelationSet`:=set(STDeclaredFrameRelation): `type/STRotationListAxis`:=[STFrameName,STFrameName,STAxisInteger,STAngle]: `type/STRotationListEvector`:=[STFrameName,STFrameName,STEvector,STAngle]: `type/STRotationListComponents`:=[STFrameName,STFrameName,{ST3VList,ST3Vector},STAngle]: `type/STRotationList`:={STRotationListAxis,STRotationListEvector,STRotationListComponents}: `type/STOrientationSpecifierSimple`:=[STFrameName,STAxisInteger,STAngle]: `type/STOrientationSpecifierAxis`:={[STFrameName,ST3Vector,STAngle],[STFrameName,ST3VList,STAngle]}: `type/STOrientationSpecifier`:={STOrientationSpecifierSimple,STOrientationSpecifierAxis}: ############################################################################## # ############################################################################## SophiaSelectEvectorComponents:=proc(evec) if (not(type([args],[STEvector]))) then ERROR(`Invalid arguments!`); fi; RETURN(evec[1]); end: ############################################################################## # ############################################################################## SophiaSelectEvectorFrame:=proc(evec) if (not(type([args],[STEvector]))) then ERROR(`Invalid arguments!`); else RETURN(evec[2]); fi; end: ############################################################################## # SophiaConstructEvector(compv,fN) # Constructs the Sophia Evector object from a list of three # components and a framename. ############################################################################## SophiaConstructEvector:=proc(compv,fN) if (not(type([args],[{ST3Vector,ST3VList},STFrameName]))) then ERROR(`Invalid arguments!`); else RETURN([convert(compv,array),fN]); fi; end: ############################################################################## # # AAAAAAA L I AAAAAAA SSSSSSS EEEEEEE SSSSSSS # A A L I A A S E S # A A L I A A S E S # AAAAAAA L I AAAAAAA SSSSSSS EEEEE SSSSSSS # A A L I A A S E S # A A L I A A S E S # A A LLLLLLL I A A SSSSSSS EEEEEEE SSSSSSS # ############################################################################## # Rmx(fN1,fN2) # Selects a rotation matrix ############################################################################## Rmx:=proc(fN1,fN2) RETURN(eval(SophiaSelectRmx(fN1,fN2))); end: ############################################################################## # # EEEEEEE X X TTTTTTT RRRRRR AAAAAAA # E X X T R R A A # E X X T R R A A # EEEEE X T RRRRRR AAAAAAA # E X X T R R A A # E X X T R R A A # EEEEEEE X X T R R A A # ############################################################################## # Functions for partitioned vectors # © 1994-1998 # Anders Lennartsson ############################################################################## SophiaGlobalAVTable:=table(): SophiaGlobalAVcTable:=table(): #################################################################### # modified angular velocity of FB in FA # # stores calculated angular velocity in table SophiaGlobalAVTable ################################################################### aV3 := proc(FA,FB) global SophiaGlobalAVTable; local b1,b2,b3,b1d,b2d,b3d,wFAFB; if type([args],[STFrameName,STFrameName]) then if ( (SophiaCheckIfRmx(FA,FB)) or (nops(SophiaSelectFindDeclaredFramePath(FA,FB))<>0) ) then b1 := Evector(1,0,0,FB); b2 := Evector(0,1,0,FB); b3 := Evector(0,0,1,FB); b1d := diffFrameTime(b1,FA); b2d := diffFrameTime(b2,FA); b3d := diffFrameTime(b3,FA); wFAFB:=Esimplify(Evector(Edot(b2d,b3),Edot(b3d,b1),Edot(b1d,b2),FB)); SophiaGlobalAVTable[FB,FA]:=Esm(-1,wFAFB); SophiaGlobalAVTable[FA,FB]:=wFAFB; RETURN(wFAFB); else ERROR(cat(`No declared path between `,FA,`and`,FB,`.`)); fi; else ERROR(`Invalid arguments!`); fi; end: #################################################################### # modified angular velocity of fB in fA # # stores calculated angular velocity in table SophiaGlobalAVcTable ################################################################### aV3c := proc(fA,fB) global SophiaGlobalAVcTable; local b1,b2,b3,b1d,b2d,b3d,wfAfB; if (not(type([args],[STFrameName,STFrameName]))) then ERROR(`Invalid arguments!`); fi; if ( (SophiaCheckIfRmx(fA,fB)) or (nops(SophiaSelectFindDeclaredFramePath(fA,fB))<>0) ) then b1 := Evector(1,0,0,fB); b2 := Evector(0,1,0,fB); b3 := Evector(0,0,1,fB); b1d := diffFrameTime(b1,fA); b2d := diffFrameTime(b2,fA); b3d := diffFrameTime(b3,fA); wfAfB:=[Esimplify(Evector(Edot(b2d,b3),Edot(b3d,b1),Edot(b1d,b2),fB))]; SophiaGlobalAVcTable[fB,fA]:=csm(-1,wfAfB); SophiaGlobalAVcTable[fA,fB]:=wfAfB; else ERROR(cat(`No declared path between `,fA,`and`,fB,`.`)); fi; end: #################################################################### # calculates time derivative of chain of Evectors relative to # specified frame ################################################################### cdft := proc(evs,fName) global SophiaGlobalAVTable; local each,evst,ltd,gtd,lfr; evst:=[]; for each in evs do ltd:=diffVectorTime(each); lfr:=SophiaSelectEvectorFrame(each); if evalb(fName<>lfr) then if (not(assigned(SophiaGlobalAVTable[fName,lfr]))) then aV3(fName,lfr); fi; gtd:=SophiaGlobalAVTable[fName,lfr] &xx each; evst:=[op(evst),ltd,gtd]; else evst:=[op(evst),ltd]; fi; od; RETURN(evst); end: #################################################################### # calculates time derivative of chain of Evectors relative to # specified frame ################################################################### cdftc := proc(evs,fName) global SophiaGlobalAVcTable; local each,evst,ltd,gtd; evst:=[]; for each in evs do ltd:=diffVectorTime(each); if evalb(fName<>(SophiaSelectEvectorFrame(each))) then if (not(assigned(SophiaGlobalAVcTable[fName,SophiaSelectEvectorFrame(each)]))) then aV3c(fName,SophiaSelectEvectorFrame(each)); fi; gtd:=ccross(SophiaGlobalAVcTable[fName,SophiaSelectEvectorFrame(each)],[each]); evst:=[op(evst),ltd,op(gtd)]; else evst:=[op(evst),ltd]; fi; od; evst; end: #################################################################### # make chain of Evectors # ################################################################### mkc := proc() [args]; end: #################################################################### # contraction of dyads with partitioned vector ################################################################### ddot:=proc(ed,ce) local ii,nc; nc:=[]; for ii in ce do nc:=[op(nc),Edot(ed,ii)]; od; RETURN(nc); end: #################################################################### # dotproduct with chain of Evectors # ################################################################### cdot := proc(c1,c2) local each1,each2,s; if type([args],[list(STEvector),list(STEvector)]) then s:=0; for each1 in c1 do for each2 in c2 do s:=s+Edot(each1,each2); od; od; RETURN(s); else ERROR(`Invalid arguments!`); fi; end: `&co`:=proc(c1,c2) cdot(c1,c2) end: #################################################################### # crossproduct with chain of Evectors # ################################################################### ccross := proc(c1,c2) local each1,each2,s; if type([args],[list(STEvector),list(STEvector)]) then s:=[]; for each1 in c1 do for each2 in c2 do s:=[op(s),Ecross(each1,each2)]; od; od; RETURN(s); else ERROR(`Invalid arguments!`); fi; end: #################################################################### # dotproduct with Kc # ################################################################### Kcdot := proc(Kc1,Kc2) local each1,each2,s,n,ii; s:=0; n:=nops(Kc1); for ii from 1 to n do s:=s+cdot(Kc1[ii],Kc2[ii]); od; s; end: #################################################################### # chain scalar multiply # ################################################################### csm := proc(s,ch) local each,nch; nch:=[]; for each in ch do nch:=[op(nch),Esm(s,each)]; od; nch; end: #################################################################### # chain express as one Evector in fN # ################################################################### cexp := proc(ch,fN) local ev,each; ev:=Evector(0,0,0,fN); for each in ch do ev:=ev &++ express(each,fN); od; RETURN(ev); end: #################################################################### # chain compact # get rid of zero vectors and move all parts in the same frame to # one Evector ################################################################### ccpt := proc(ch,fN) local each,nch,evtable,frameset,eachf; nch:=[]; frameset:={}; for each in ch do eachf:=SophiaSelectEvectorFrame(each); if not(member(eachf,frameset)) then frameset:=frameset union {SophiaSelectEvectorFrame(each)}; evtable[eachf]:=each; else evtable[eachf]:=evtable[eachf] &++ each; fi; od; for eachf in frameset do each:=evtable[eachf]; if not(Ec(each,1)=0 and Ec(each,2)=0 and Ec(each,3)=0) then nch:=[op(nch),each]; fi; od; if nch=[] then nch:=[Evector(0,0,0,eachf)]; fi; RETURN(nch); end: #################################################################### # Kchain scalar multiply # ################################################################### Kcsm := proc(s,Kc1) local Kc,n,ii; Kc:=[]; n:=nops(Kc1); for ii from 1 to n do Kc:=[op(Kc),csm(s,ccpt([op(Kc1[ii])]))]; od; Kc; end: #################################################################### # chain add ################################################################### cadd := proc() local ii; RETURN([seq(op(args[ii]),ii=1..nargs)]); end: #################################################################### # Kchain add # ################################################################### Kcadd := proc(Kc1,Kc2) local Kc,n,ii; Kc:=[]; n:=nops(Kc1); for ii from 1 to n do Kc:=[op(Kc),ccpt([op(Kc1[ii]),op(Kc2[ii])])]; od; Kc; end: #################################################################### # Kchain subtract # ################################################################### Kcsub := proc(Kc1,Kc2) local ec,Kc,n,ii; Kc:=[]; n:=nops(Kc1); for ii from 1 to n do Kc:=[op(Kc),ccpt([op(Kc1[ii]),op(csm(-1,Kc2[ii]))])]; od; Kc; end: #################################################################### # chain tangents # ################################################################### Kctau:=proc(cv,var,dimen) local ii,jj,n,imx,zs,ctvt,vvt,ctv,eachc,eachev; ctv:=[]; if type([args],[list(list(STEvector)),STuaname,integer]) then zs:={seq(cat(var,ii)=0,ii=1..dimen)}; ctvt:=map(ccpt,subs(zs,cv)); print(`Explicitly time dependent part:`,ctvt); imx:=array(1..dimen,1..dimen,identity); for ii from 1 to dimen do zs:={seq(cat(var,jj)=imx[ii,jj],jj=1..dimen)}; ctv:=[op(ctv),Kcsub(subs(zs,cv),ctvt)]; od; elif type([args],[list(list(STEvector)),list(STuaname)]) then n:=nops(var); zs:={seq(var[ii]=0,ii=1..n)}; ctvt:=map(ccpt,subs(zs,cv)); print(`Explicitly time dependent part:`,ctvt); imx:=array(1..n,1..n,identity); for ii from 1 to n do zs:={seq(var[jj]=imx[ii,jj],jj=1..n)}; ctv:=[op(ctv),Kcsub(subs(zs,cv),ctvt)]; od; else ERROR(`Invalid arguments!`); fi; RETURN(ctv); end: #################################################################### # chain kane equations # ################################################################### ckane:=proc(taus,fa,pt) local keqs,nn,ii; keqs:=[]; nn:=nops(taus); for ii from 1 to nn do keqs:=[op(keqs),Kcdot(taus[ii],fa)-Kcdot(taus[ii],pt)]; od; keqs; end: #################################################################### # SophiaClearGlobals() # Clear global variables. #################################################################### SophiaClearGlobals:=proc() global SophiaGlobalAVTable, SophiaGlobalAVcTable, SophiaGlobalRmxTable, SophiaGlobalDeclaredFrameRelations, SophiaGlobalRmxTableApprox, SophiaGlobalRmxTableExact; SophiaGlobalRmxTable:=table(); SophiaGlobalRmxTableApprox:='SophiaGlobalRmxTableApprox'; SophiaGlobalRmxTableExact:='SophiaGlobalRmxTableExact'; SophiaGlobalAVTable:=table(); SophiaGlobalAVcTable:=table(); SophiaGlobalDeclaredFrameRelations:={}; print(`Global variables cleared!`); end: #################################################################### # outer product for partitioned vectors coutp := proc(c1,c2) local each1,each2,s; if type([args],[list(STEvector),list(STEvector)]) then s:=[]; for each1 in c1 do for each2 in c2 do s:=[op(s),Eoutp(each1,each2)]; od; od; RETURN(s); else ERROR(`Invalid arguments!`); fi; end: