In questa sezione riportiamo il sorgente del pacchetto NImmersion.m
per la parte riguardante il codice delle funzioni, per avere un'utile
refferenza su carta del codice del programma. La versione completa si può
comunque trovare sul disco in allegato.
Download Notebook NImmersion.ma
Off[General::spell1]
HaccaK1=.
dg=.
opts=.
opt=.
gopts=.
dxsur=.
dtsur=.
NHaccaK1=.
ires=.
sum=.
res=.
ddt=.
ddtt=.
tpoints=.
ntpoints=.
nxpoints=.
npoints=.
th0=.
th1=.
th2=.
th0t=.
th1t=.
th2t=.
th1x=.
th0x=.
tmp=.
temp=.
tmp1=.
tmp2=.
tmpr=.
tmpv=.
tmpg=.
on=.
og=.
RVector=.
VVector=.
tsur=.
sur=.
SamplePoints$x=.
SamplePoints$t=.
ConstrainAngle1=.
ConstrainAngle2=.
xpoints=.
On[General::spell1]
(*
:[font = input; initialization; preserveAspect]
*)
Options[ND]^={WorkingPrecision->$MachinePrecision};
ND[data_,h_,opts___]:=
Block[{l=Length[data],r,i,p},
p=WorkingPrecision /. {opts} /. Options[ND];
r=If[ l>2,
N[(Drop[data,2]-Drop[data,-2])/(2*h),p]
(* else *),
{}
];
(* endif *)
(* Agli estremi la derivata \`e calcolata *)
(* In maniera tradizionale *)
AppendTo[r,N[(data[[l]]-data[[l-1]])/h,p]];
PrependTo[r,N[(data[[2]]-data[[1]])/h,p]]
/; ( NumberQ[h] && VectorQ[data,NumberQ] && l > 1)
];
(*
:[font = input; initialization; preserveAspect]
*)
ND[data_,opts___]:=
Block[{l=Length[data],r,rh,s,p},
p=WorkingPrecision /. {opts} /. Options[ND];
s=If[ l>2,
{rh,r}=Transpose[data];
N[(Drop[r,2]-Drop[r,-2])/
(Drop[rh,2]-Drop[rh,-2]),p]
(* else *),
{}
];
(* endif *)
(* Agli estremi la derivata \`e calcolata *)
(* In maniera tradizionale *)
AppendTo[s,N[(data[[l]][[2]]-data[[l-1]][[2]])/
(data[[l]][[1]]-data[[l-1]][[1]]),p]];
PrependTo[s,N[(data[[2]][[2]]-data[[1]][[2]])/
(data[[2]][[1]]-data[[1]][[1]]),p]];
Transpose[{rh,s}]
/; ( MatrixQ[data,NumberQ] && l > 1 &&
Length[ data[[1]] ] == 2 )
];
(*
:[font = input; initialization; preserveAspect]
*)
Options[MakeData]^={WorkingPrecision->$MachinePrecision,
SamplePoints->50,
Compiled->True};
MakeData::ppts="SamplePoints options must be greater than 1";
MakeData::prcsm="WorkingPrecision must be greater than 0";
MakeData[f_,x_Symbol,{xmin_,xmax_},opts___]:=
Block[{fc,y,p,n},
p=WorkingPrecision /. {opts} /. Options[MakeData];
n=SamplePoints /. {opts} /. Options[MakeData];
If[!NumberQ[n] || n<2,
Message[MakeData::ppts];Return[$Failed]];
If[!NumberQ[p] || p<1,
Message[MakeData::prcsm];Return[$Failed]];
If[p<=$MachinePrecision &&
(Compiled /. {opts} /. Options[MakeData]),
(*Then*) fc=Compile[x, f],
(*Else*) fc=Function[x,f]];
Table[{y,fc[y]},
{y,N[xmin,p],N[xmax,p],
N[(xmax-xmin)/(n-1),p]}]
];
(*
:[font = input; initialization; preserveAspect]
*)
HaccaK1[t_Symbol,cet_,cct_,k1_]:=(D[cct,t]/(k1*cet));
(*
:[font = input; initialization; preserveAspect]
*)
Options[NHaccaK1]^={WorkingPrecision->$MachinePrecision};
NHaccaK1[dcct_,dcet_,k1_,opts___]:=
Block[{p,temp},
p=WorkingPrecision /. {opts} /. Options[NHaccaK1];
Transpose[{temp[[1]],
N[(Transpose[ND[dcct,opts]][[2]]/
(k1*temp[[2]])),p]}
] /; (Length[dcet]==Length[dcct] &&
(temp=Transpose[dcet])[[1]]==
Transpose[dcct][[1]] &&
AtomQ[k1])];
(*
:[font = input; initialization; preserveAspect]
*)
(*Definizione per i casi particolari *)
InverseSin[Sin[h_],ep_]:=h /; ep==1;
InverseSin[Sinh[h_],ep_]:=h /; ep==-1;
InverseSin[Cos[h_],ep_]:=Pi/2-h /; ep==1;
InverseSin[-Sin[h_],ep_]:=-h /; ep==1;
InverseSin[-Sinh[h_],ep_]:=-h /; ep==-1;
InverseSin[-Cos[h_],ep_]:=-Pi/2+h /; ep==1;
(*
:[font = input; initialization; preserveAspect]
*)
Options[OptimizeDelta]^={
StartSamplePoints->10,
WorkingPrecision->$MachinePrecision,
Compiled->True,
MaxRecursion->6,
Subdivisions->2
};
(*
:[font = input; initialization; preserveAspect]
*)
OptimizeDelta::fewp="StartSamplePoints must be greater than 1";
OptimizeDelta[f_,x_Symbol,{x0_,x1_},maxinc_Real,opts___]:=
Block[{val1,val2,val3,delta,l,m,i,st,en,div,fc,startn,sudd},
p=WorkingPrecision /. {opts} /. Options[OptimizeDelta];
If[p<=$MachinePrecision &&
(Compiled /. {opts} /. Options[OptimizeDelta]),
(*Then*) fc=Compile[x, f],
(*Else*) fc=Function[x,f]];
If[(startn=(StartSamplePoints /.
{opts} /. Options[OptimizeDelta]))<=2,
(*Then*) Message[OptimizeDelta::fewp];Return[$Failed]];
sudd=Subdivisions /. {opts} /. Options[OptimizeDelta];
--startn;
delta=N[(x1-x0)/(startn),p];
val1=Table[N[fc[x],p] , {x, x0, x1, delta}];
val2=ND[val1,delta,opts];
val3=Table[x , {x, x0, x1, delta}];
l=Length[val2];
i=Indeterminate;
While[ m=-Infinity;
Do[ If[m<(val2[[j]]*delta), m=val2[[j]]*delta;i=j], {j,l}];
m>maxinc,
(*Do*)
Which[i==1,
st=1;en=2,
i==l,
st=i-1;en=i,
i>1 && i<l,
st=i-1;en=i+1,
True,
st=Max[i-1,1];en=Min[l,i+1]
];
delta/=sudd;
val1=Table[N[fc[x],p] , {x, val3[[st]],val3[[en]], delta}];
val2=ND[val1,delta,opts];
l=Length[val2];
val3=Table[N[x,p] , {x, val3[[st]],val3[[en]], delta}]
];
startn=Floor[N[(x1-x0)/delta,p]]+1;
{startn,N[(x1-x0)/(startn),p]}
] /; NumberQ[x0] && NumberQ[x1] && x0<x1;
(*
:[font = input; initialization; preserveAspect]
*)
Options[ImmersionDecomponible]^={
DefaultSign->1,
MetricESign->1,
WorkingPrecision->$MachinePrecision,
Kappa1->1
};
(*
:[font = input; initialization; preserveAspect]
*)
ImmersionDecomponible[cet_, ccx_, cct_,
{t_Symbol, x_Symbol},{a_,b_},sur_,fbb_,opts___]:=
Block[{fa,fb,ifb,dfb,h,res,ires,res1,res2,ep,sgn,k1},
ep=MetricESign /. {opts} /. Options[ImmersionDecomponible];
sgn=DefaultSign /. {opts} /. Options[ImmersionDecomponible];
k1=Kappa1 /. {opts} /. Options[ImmersionDecomponible];
(*Calcola l'angolo \Theta^1*)
ifb=k1*Integrate[ccx,x];
fb=ifb - ((ifb) /. x->b)+fbb;
(*Calcola l'angolo \Theta^0 se possibile*)
fa=InverseSin[HaccaK1[t,cet,cct,k1],ep];
(*Calcola le componenti del risultato*)
res=If[ MatchQ[fa,InverseSin[__]],(*Controlla se \`e un caso particolare*)
(*Then*) If[ep==-1,1,sgn]/Abs[k1] *
Sqrt[ (k1*cet)^2 - ep*D[cct,t]^2 ],
(*Else*) If[ep==1,cet*Cos[fa],cet*Cosh[fa]]
(*EndIf*) ];
ires=Integrate[res,t];
res1=cct/k1*Sin[fb];
res2=cct/k1*Cos[fb];
(*Vettore del risultato*)
{ires-((ires) /. t->a) + sur[[1]],
res1-((res1) /. {t->a,x->b}) + sur[[2]],
res2-((res2) /. {t->a,x->b}) + sur[[3]]}
] /; (FreeQ[ccx,t] && FreeQ[{cct,cce},x] &&
AtomQ[a] && AtomQ[b] && AtomQ[fbb] &&
VectorQ[sur] && Length[sur]==3);
(*\scriptsize\end{verbatim} *)
(*
:[font = input; initialization; preserveAspect]
*)
ListInt2[cl_,k_] :=
Block[{i,j,l,tmp,sum,x,m=k/2,m1=k/2+1},
sum=Table[0,{i,(l=Length[cl])}];
If[Length[cl] == k,
Do[sum[[j]]=IntInt[cl,1,j],{j,2,k}],
(* else *)
Do[sum[[j]]=IntInt[cl,1,j],{j,2,m1}];
Do[sum[[m+i]] = sum[[m+i-1]]+
IntInt[Take[cl,{i,i+k-1}],m,m1],
{i,2,l-k}
];
tmp = Take[cl,-k];
Do[sum[[l-k+j]]=sum[[l-k+m]]+
IntInt[tmp,m,j],{j,m1,k}];
];
Table[{cl[[i]][[1]],sum[[i]]},{i,l}]
];
(*
:[font = input; initialization; preserveAspect]
*)
IntInt[data_,i_,j_] :=
Block[{cl, x, k, sum, tmp},
cl = CoefficientList[InterpolatingPolynomial[data,x],x];
x = data[[j,1]];
tmp = 0;
Do[tmp = (tmp + cl[[k]]/k) x,{k,Length[cl],1,-1}];
sum = tmp;
x = data[[i,1]];
tmp = 0;
Do[tmp = (tmp + cl[[k]]/k) x,{k,Length[cl],1,-1}];
sum - tmp
];
(*
:[font = input; initialization; preserveAspect]
*)
ListIntegrate2[cl_?VectorQ, h_, k_Integer:4] :=
Block[{ans},
ans /; (ans = ListIntegrate0[h, cl, k]) =!= $Failed
] /; k > 0;
ListIntegrate2[cl_?MatrixQ, k_Integer:4] :=
Block[{ans, dim = Length[Dimensions[cl]]},
ans /; ((dim == 2) &&
(Length[cl[[1]]] == 2) &&
((ans = ListIntegrate1[cl, k]) =!= $Failed))
] /; k > 0;
ListIntegrate::degerr = "Number of data points (`1`) is insufficient for the
degree of interpolant requested (`2`).";
ListIntegrate0[h_, cl_, k_] :=
Block[{ke = k + If[IntegerQ[k/2], 0, 1], data,l,sum,y},
If[Length[cl] < ke,
Message[ListIntegrate::degerr,Length[cl],ke-1];
Return[$Failed]
];
If[ke == 2,
l = Length[cl];
y = Drop[cl,1]+Drop[cl,-1];
sum=Table[0,{i,l}];
Do[sum[[i+1]]=sum[[i]]+h y[[i]]/2,{i,l-1}];
Table[{h*i,sum[[i]]},{i,l}],
(* else *)
data = Array[h #&,{Length[cl]}];
data = Transpose[{data,cl}];
ListInt2[data,ke]
]
];
ListIntegrate1[cl_, k_] :=
Block[{ke = k + If[IntegerQ[k/2], 0, 1], x, y,l,sum},
If[Length[cl] < ke,
Message[ListIntegrate::degerr,Length[cl],ke-1];
Return[$Failed]
];
If[ke == 2,
y = Transpose[cl];
l = Length[cl];
x = Drop[y[[1]],1]-Drop[y[[1]],-1];
y = Drop[y[[2]],1]+Drop[y[[2]],-1];
sum=Table[0,{i,l}];
Do[sum[[i+1]]=sum[[i]]+x[[i]]y[[i]]/2,{i,l-1}];
Table[{cl[[i]][[1]],sum[[i]]},{i,l}],
(* else *)
ListInt2[cl,ke]
]
];
(*
:[font = input; initialization; preserveAspect]
*)
EFGComponents[sur:{_,_,_},{t_Symbol,x_Symbol},1]:=
Block[{dx,dt},dt=D[sur,t];dx=D[sur,x];
{dt.dt,dt.dx,dx.dx} // Simplify
]
(*
:[font = input; initialization; preserveAspect]
*)
EFGComponents[sur:{_,_,_},{t_Symbol,x_Symbol},-1]:=
Block[{dx2,dx,dt2,dt},dt=D[sur,t];dx=D[sur,x];
dt2={-dt[[1]],dt[[2]],dt[[3]]};
dx2={-dx[[1]],dx[[2]],dx[[3]]};
{dt.dt2,dt2.dx,dx2.dx} // Simplify
]
(*
:[font = input; initialization; preserveAspect]
*)
EpQ=(NumberQ[#] && (#==1 || #==-1))&;
(*
:[font = input; initialization; preserveAspect]
*)
Kurvature[1,g_,{t_Symbol,x_Symbol},ep_?EpQ]:=
Simplify[-1/g *
ep*D[g,{t,2}]];
(*
:[font = input; initialization; preserveAspect]
*)
Kurvature[e_,g_,{x_Symbol,t_Symbol},ep_?EpQ]:=
Simplify[-1/(e*g) *
ep*D[Evaluate[D[g,t]/e],t] +
D[Evaluate[D[e,x]/g],x]];
(*
:[font = input; initialization; preserveAspect]
*)
ColorStretched[v_?NumberQ,s_?NumberQ]=
Hue[ArcTan[v s]/Pi+0.5];
ColorStretched[v_/;(v==Infinity || v==-Infinity ||
v==ComplexInfinity), s_]=Hue[0];
ColorStretched[v_, s_/;(s==Infinity || s==-Infinity ||
s==ComplexInfinity)]=Hue[0];
(*
:[font = input; initialization; preserveAspect]
*)
PlotLegend[vmin_,vmax_,off_,s_,p_:100]:=Block[{nil},
DensityPlot[ArcTan[(v-off) s]/Pi+0.5,{v,vmin,vmax},{nil,-2,2},
Axes->{True,False},
AspectRatio->1/4,
RotateLabel->False,
FrameTicks->{Automatic,None},
AxesLabel->{"Value",""},
PlotPoints->{p,2},
Mesh->False,
AspectRatio->1/4,
PlotLabel->"Curvatura",
ColorFunction->Hue
]
];
(*
:[font = input; initialization; preserveAspect]
*)
InverseCosSin[a_?VectorQ,b_?VectorQ,
st_?NumberQ,wp:{_,_}:{10^-10,$MachinePrecision}]:=
Block[{vecchio,i,l,ris,rot,ang,z,w,p},(
p=wp[[2]];w=N[wp[[1]],p];
l=Length[a];
ang=N[Pi/2,p];
{ris[0],rot[0]}=ConstrainAngle2[st,p];
Do[ z=N[a[[i]]+I b[[i]], p];
If[Chop[z,w]==0,
(*Then*) ris[i]=ris[i-1];rot[i]=rot[i-1],
(*Else*) ris[i]=N[Arg[ z ],p];
Which[ris[i-1]<(-ang) && ris[i]>ang,
rot[i]=rot[i-1]-1,
ris[i-1]>ang && ris[i]<(-ang),
rot[i]=rot[i-1]+1,
True,
rot[i]=rot[i-1]
]
(*EndIf*)],
{i,1,l}];
{Table[ris[i],{i,l}],Table[rot[i],{i,l}]}
) /; (l=Length[a])==Length[b] ];
(*
:[font = input; initialization; preserveAspect]
*)
ConstrainAngle2[a_?NumberQ,p_:$MachinePrecision]:=
Block[{rot,c},
rot=N[Quotient[a,2 Pi],p];
If[rot<0,++rot];
c=N[a-rot*2 Pi,p];
If[c>N[Pi,p] || c<N[-Pi,p], c=N[c- Sign[c] 2 Pi,p];rot=rot-Sign[c] ];
{N[a-rot*2 Pi,p],rot}
];
SetAttributes[ConstrainAngle2,Listable];
(*
:[font = input; initialization; preserveAspect]
*)
ConstrainAngle1[a_?NumberQ,p_:$MachinePrecision]:=
Block[{rot,c},
rot=N[Quotient[a,2 Pi],p];
If[rot<0,++rot];
c=N[a-rot*2 Pi,p];
If[c>N[Pi,p] || c<N[-Pi,p], c=N[c- Sign[c] 2 Pi,p];rot=rot-Sign[c] ];
N[a-rot*2 Pi,p]
];
SetAttributes[ConstrainAngle1,Listable];
(*
:[font = input; initialization; preserveAspect]
*)
ParametricListPlot3D::wrdim=
"All the matrix must have the same dimension";
(*
;[s]
1:0,0;77,-1;
1:1,0,0 ,Courier,0,10,0,0,0;
:[font = input; initialization; preserveAspect]
*)
ParametricListPlot3D[
sur1_?MatrixQ,
sur2_?MatrixQ,
sur3_?MatrixQ,
color_?MatrixQ,opts___]:=
Block[{ntpoints,nxpoints,a,b,c,d},
If[!(Dimensions[sur1]==Dimensions[sur2]==Dimensions[sur3]),
Message[ParametricListPlot3D::wrdim];Return[$Failed];];
{nxpoints,ntpoints}=Dimensions[sur1];
Show[Graphics3D[Flatten[Table[
a=Sequence[n-1,j-1];
b=Sequence[n-1,j];
c=Sequence[n,j];
d=Sequence[n,j-1];
{color[[a]],Polygon[
{{sur1[[a]],sur2[[a]],sur3[[a]]},
{sur1[[b]],sur2[[b]],sur3[[b]]},
{sur1[[c]],sur2[[c]],sur3[[c]]},
{sur1[[d]],sur2[[d]],sur3[[d]]}}
]},
{n,2,nxpoints},{j,2,ntpoints}] ]], Lighting->False,opts ]
];
(*
;[s]
1:0,0;649,-1;
1:1,0,0 ,Courier,0,10,0,0,0;
:[font = input; initialization; preserveAspect]
*)
ParametricListPlot3D[
sur1_?MatrixQ,
sur2_?MatrixQ,
sur3_?MatrixQ,opts___]:=
Block[{ntpoints,nxpoints,a,b,c,d},
If[!(Dimensions[sur1]==Dimensions[sur2]==Dimensions[sur3]),
Message[ParametricListPlot3D::wrdim];Return[$Failed];];
{nxpoints,ntpoints}=Dimensions[sur1];
Show[Graphics3D[Flatten[Table[
a=Sequence[n-1,j-1];
b=Sequence[n-1,j];
c=Sequence[n,j];
d=Sequence[n,j-1];
{Polygon[
{{sur1[[a]],sur2[[a]],sur3[[a]]},
{sur1[[b]],sur2[[b]],sur3[[b]]},
{sur1[[c]],sur2[[c]],sur3[[c]]},
{sur1[[d]],sur2[[d]],sur3[[d]]}}
]},
{n,2,nxpoints},{j,2,ntpoints}] ]],opts ]
];
(*
:[font = input; initialization; preserveAspect]
*)
Options[ThetaInit]^={
SamplePoints$x->50,
SamplePoints$t->50,
Precision->$MachinePrecision,
WarningLimit->0.001
};
(*
:[font = input; initialization; preserveAspect]
*)
ThetaInit::bopt="Bad options, Samplepoints must be >1 and
Precision must be > 0,WarningLimit must be a number.";
ThetaInit::crit="Warning, critical point `2` in `1`";
ThetaInit::errr="Calculus of `1` is broken";
(*
:[font = input; initialization; preserveAspect]
*)
ThetaInit[1,g_,{t_Symbol,tmin_?NumberQ,tmax_?NumberQ},
{x_Symbol,xmin_?NumberQ,xmax_?NumberQ},
opts___]:=
Block[{msg,m,w,x0,p,tpoints,npoints,fa,fb,fs,hk1,
i,ddt,crit,ddtt},
p=Precision /. {opts} /. Options[ThetaInit];
w=WarningLimit /. {opts} /. Options[ThetaInit];
npoints=SamplePoints /. {opts} /. Options[ThetaInit];
If[ !VectorQ[{p,npoints,w},NumberQ] ||
p<1 || npoints<2 ,
(*Then*) Message[ThetaInit::bopt];Return[$Failed]];
x0=N[xmin,p];
ddt=D[ Evaluate[(g /. x->x0)] ,t];
If[FreeQ[ddt,Derivative],
(*Then*) ddt=MakeData[ddt,t,{tmin,tmax},opts];,
(*Else*) ddt=ND[MakeData[(g /. x->x0),
t,{tmin,tmax},opts],opts];
(*EndIf*)];
If[!MatrixQ[ddt,NumberQ],
msg=StringForm["D[`1`,`2`]",g,t];
Message[ThetaInit::errr,msg];Return[$Failed]];
crit=Select[ddt,(Abs[#[[2]]]<Abs[w])&];
If[Length[crit]>0,
Message[ThetaInit::crit,Short[ Transpose[crit][[1]] ],
"-Possible numerical instability- near t"]
];
ddtt=Transpose[ddt];
If[ (m=Max[Abs[ ddtt[[2]] ] ])>1,
(*Then*)
hk1=Transpose[{ddtt[[1]],N[ddtt[[2]]/m,p]}],
(*Else*)
hk1=ddt;m=1;
];
crit=Select[hk1,(Abs[#[[2]]]>Abs[1-Abs[w]])&];
If[Length[crit]>0,
Message[ThetaInit::crit,Short[ Transpose[crit][[1]] ],
"-Possible inderivability- near t"]
];
hk1=Transpose[hk1];
tpoints=hk1[[1]];
fa=N[ArcSin[hk1[[2]]],p];
fb=N[m*x0,p];
fs=N[Pi/2,p];
If[!VectorQ[fa,NumberQ] || !NumberQ[fb] ||
!NumberQ[fs] ,
Message[ThetaInit::errr,"initial angles"];Return[$Failed]];
{x0,tpoints,{
fa,
Table[fb,{i,npoints}],
Table[fs,{i,npoints}]
}
}
] /; xmax > xmin && tmax > tmin;
(*
:[font = input; initialization; preserveAspect]
*)
ThetaInit[-1,g_,{t_Symbol,tmin_?NumberQ,tmax_?NumberQ},
{x_Symbol,xmin_?NumberQ,xmax_?NumberQ},
opts___]:=
Block[{msg,m,w,x0,p,npoints,fa,fb,fs,ddt,crit,ddtt},
p=Precision /. {opts} /. Options[ThetaInit];
w=WarningLimit /. {opts} /. Options[ThetaInit];
npoints=SamplePoints /. {opts} /. Options[ThetaInit];
If[ !VectorQ[{p,npoints,w},NumberQ] ||
p<1 || npoints<2 ,
(*Then*) Message[ThetaInit::bopt];Return[$Failed]];
x0=N[xmin,p];
ddt=D[Evaluate[(g /. x->x0)],t];
If[FreeQ[ddt,Derivative],
(*Then*) ddt=MakeData[ddt,t,{tmin,tmax},opts],
(*Else*) ddt=ND[MakeData[(g /. x->x0)
,t,{tmin,tmax},opts],opts];
(*EndIf*)];
If[!MatrixQ[ddt,NumberQ],
msg=StringForm["D[`1`,`2`]",g,t];
Message[ThetaInit::errr,msg];Return[$Failed]];
crit=Select[ddt,(Abs[#[[2]]]<Abs[w])&];
If[Length[crit]>0,
Message[ThetaInit::crit,Short[ Transpose[crit][[1]] ],
"-Possible numerical instability- for t"]
];
ddtt=Transpose[ddt];
tpoints=ddtt[[1]];
fa=N[ArcSinh[ddtt[[2]]],p];
fb=N[x0,p];
fs=N[Pi/2,p];
If[!VectorQ[fa,NumberQ] || !NumberQ[fb] ||
!NumberQ[fs] ,
Message[ThetaInit::errr,"initial angles"];Return[$Failed]];
{x0,tpoints,{
fa,
Table[fb,{i,npoints}],
Table[fs,{i,npoints}]
}
}
] /; xmax > xmin && tmax > tmin;
(*
:[font = input; initialization; preserveAspect]
*)
thetax[th0_,th2_,th1t_,th2t_,gt_,g_,ep_]:=
Block[{thh,th0x,th1x},
thh= g (th2t - th1t Cos[th0]);
th0x= -gt Cos[th2] + thh Sin[th2];
th1x= Csc[th0] (thh Cos[th2] + gt Sin[th2]);
{th0x,th1x}
] /;ep==1;
(*
:[font = input; initialization; preserveAspect]
*)
thetax[th0_,th2_,th1t_,th2t_,gt_,g_,ep_]:=
Block[{thh,th0x,th1x},
thh= g (th2t + th1t Cosh[th0]);
th0x= gt Cos[th2] - thh Sin[th2];
th1x= Csch[th0] (thh Cos[th2] + gt Sin[th2]);
{th0x,th1x}
] /;ep==-1;
(*
:[font = input; initialization; preserveAspect]
*)
RVector[f1_,f2_,ep_]:={
Cos[f1],
Sin[f1]*Sin[f2],
Sin[f1]*Cos[f2]
} /;ep==1; (*Modulo 1*)
VVector[f1_,f2_,fs_,ep_]:={
Cos[fs] Sin[f1],
-Cos[fs] Cos[f1] Sin[f2] + Sin[fs] Cos[f2], (*Modulo 1*)
-Cos[fs] Cos[f1] Cos[f2] - Sin[fs] Sin[f2]
} /;ep==1;
RVector[f1_,f2_,ep_]:={
Cosh[f1],
Sinh[f1]*Sin[f2], (*Modulo -1*)
Sinh[f1]*Cos[f2]
} /;ep==-1;
VVector[f1_,f2_,fs_,ep_]:={
Cos[fs] Sinh[f1],
Cos[fs] Cosh[f1] Sin[f2] + Sin[fs] Cos[f2], (*Modulo 1*)
Cos[fs] Cosh[f1] Cos[f2] - Sin[fs] Sin[f2]
} /;ep==-1;
SetAttributes[RVector,Listable];
SetAttributes[VVector,Listable];
(*
:[font = input; initialization; preserveAspect]
*)
Options[ImmersionOf]^={
SamplePoints$t->50,
SamplePoints$x->100,
Precision->$MachinePrecision,
WarningLimit->0.001
};
(*
:[font = input; initialization; preserveAspect]
*)
ImmersionOf::bopt="Bad options, Samplepoints must be >1 and
Precision must be > 0,WarningLimit must be a number,
SamplePoints$t>5, SamplePoints$x>2";
ImmersionOf::crit="Warning, critical point `2` in `1`";
ImmersionOf::errr="Calculus of `1` is broken";
(*
:[font = input; initialization; preserveAspect]
*)
ImmersionOf[g_,tin:{t_Symbol,tmin_?NumberQ,tmax_?NumberQ},
xin:{x_Symbol,xmin_?NumberQ,xmax_?NumberQ},
{th0_Symbol,th1_Symbol,th2_Symbol},
ng_Symbol,{tdelta_Symbol, xdelta_Symbol},
ep_?EpQ,opt___]:=
Block[{j,n,i,on,iniz,tpoints,y,x0,tmp,ntpoints,p,w,
oldth0t,oldth1t,oldth2t,ogt,xpoints,nxpoints,
thx,tmp1,tmp2,crit,opts},
p=Precision /. {opt} /. Options[ImmersionOf];
w=WarningLimit /. {opt} /. Options[ImmersionOf];
ntpoints=SamplePoints$t /. {opt} /. Options[ImmersionOf];
nxpoints=SamplePoints$x /. {opt} /. Options[ImmersionOf];
If[ !VectorQ[{p,ntpoints,nxpoints,w},NumberQ] ||
p<1 || ntpoints<6 || nxpoints<2,
(*Then*) Message[ImmersionOf::bopt];Return[$Failed]];
opts=Sequence @@ Join[{SamplePoints->ntpoints},{opt}];
If[(iniz=ThetaInit[ep,g,tin,xin,opts])
==$Failed,
(*Then*) Return[$Failed]];
x0=iniz[[1]];
tpoints=iniz[[2]];
tdelta=N[(tpoints[[-1]]-tpoints[[1]])/(ntpoints-1),p];
xdelta=N[(xmax-x0)/(nxpoints-1),p];
xpoints=Table[N[y,p],{y,xmin,xmax,xdelta}];
th0[1]=iniz[[3]][[1]];
th1[1]=iniz[[3]][[2]];
th2[1]=iniz[[3]][[3]];
thx[j_,n_]:=thetax[
th0[n][[j]],
th2[n][[j]],
oldth1t[[j]],oldth2t[[j]],
ogt[[j]],ng[n][[j]],ep];
oldth0t=ND[th0[1],tdelta,opts];
oldth1t=ND[th1[1],tdelta,opts];
oldth2t=ND[th2[1],tdelta,opts];
Do[
on=n-1;
ng[on]=Table[N[(g /. {x->xpoints[[on]],t->y}), p],
{y,tmin,tmax,tdelta}];
ogt=ND[ng[on],tdelta,opts];
(*kurv[on]=-ep*ND[ogt,tdelta,opts]/ng[on];*)
tmp=Table[N[thx[j,on],p], {j,2,ntpoints-1}];
tmp=Transpose[tmp];
th0[n]=Table[(th0[on][[j+1]]+th0[on][[j-1]])/2,
{j,2,ntpoints-1}] + tmp[[1]]*xdelta;
th1[n]=Table[(th1[on][[j+1]]+th1[on][[j-1]])/2,
{j,2,ntpoints-1}] + tmp[[2]]*xdelta;
tmp=(InterpolatingPolynomial[th0[n][[{1,2,3,4}]],y] /. y->0);
tmp1=N[thx[1,on],p] xdelta + {th0[on][[1]],th1[on][[1]]};
PrependTo[th0[n],N[tmp,p]];
tmp=(InterpolatingPolynomial[th1[n][[{1,2,3,4}]],y] /. y->0);
PrependTo[th1[n],N[tmp,p]];
tmp=(InterpolatingPolynomial[th0[n][[{-1,-2,-3,-4}]],y] /. y->0);
tmp1=N[thx[-1,on],p] xdelta + {th0[on][[-1]],th1[on][[-1]]};
AppendTo[th0[n],N[tmp,p]];
crit=Select[th0[n],(Abs[#]<Abs[w])&];
If[Length[crit]>0,
Message[ImmersionOf::crit,xpoints[[n]],
"-Possible numerical instability- near x"]
];
tmp=(InterpolatingPolynomial[th1[n][[{-1,-2,-3,-4}]],y] /. y->0);
AppendTo[th1[n],N[tmp,p]];
(*Dialog[]*)
oldth0t=ND[th0[n],tdelta,opts];
oldth1t=ND[th1[n],tdelta,opts];
tmp=InverseCosSin[
oldth1t If[ep==-1,
(*Then*) Sinh[th0[n]],
(*Else*) -Sin[th0[n] ]
(*EndIf*) ],
-oldth0t,th2[on][[1]],{w,p}];
th2[n]=N[tmp[[1]]+2 Pi tmp[[2]],p];
oldth2t=ND[th2[n],tdelta,opts];,
{n,2,nxpoints}
];
on=nxpoints;
ng[on]=Table[N[(g /. {x->xpoints[[on]],t->y}), p],
{y,tmin,tmax,tdelta}];
(*ogt=ND[ng[on],tdelta,opts];*)
(*kurv[on]=-ep*ND[ogt,tdelta,opts]/ng[on];*)
] /; xmax > xmin && tmax > tmin;
(*
:[font = input; initialization; preserveAspect]
*)
ImmersedSurface[{th0_Symbol,th1_Symbol,th2_Symbol},
ng_Symbol,{tdelta_?NumberQ,xdelta_?NumberQ},
ep_?EpQ,k_Integer,sur_Symbol,opts___]:=
Block[{ntpoints,nxpoints,n,j,i,tmp,tmp2,tmpr,tmpv,
tmpg,tsur},
p=Precision /. {opts} /. Options[ImmersionOf];
ntpoints=SamplePoints$t /. {opts} /. Options[ImmersionOf];
nxpoints=SamplePoints$x /. {opts} /. Options[ImmersionOf];
tmpr=Transpose[RVector[th0[1],th1[1],ep]];
Do[(tmp[i]=Transpose[
ListIntegrate2[tmpr[[i]],tdelta,k]][[2]]),
{i,3}];
Do[
tmp2=Table[{th0[n][[j]],th1[n][[j]],th2[n][[j]]},
{n,nxpoints}];
tmp2=Transpose[tmp2];
tmpv=Transpose[VVector[Sequence @@ tmp2, ep ]];
tmpg=Table[ng[n][[j]],{n,nxpoints}];
Do[tsur[i,j]=
Transpose[ListIntegrate2[tmpg*tmpv[[i]], xdelta,k]+
tmp[i][[j]]][[2]],
{i,3}];,
{j,1,ntpoints}];
Do[sur[i]=Table[tsur[i,j][[n]],
{n,nxpoints},{j,ntpoints}],{i,3}];
];
(*
:[font = input; initialization; preserveAspect]
*)
PlotImmersionError[
sur1_?MatrixQ,
sur2_?MatrixQ,
sur3_?MatrixQ,
ng_Symbol,ep_?EpQ,{tdelta_?NumberQ,xdelta_?NumberQ},
opts___]:=
Block[{sur,ntpoints,nxpoints,a,gg2,b,c,d,mm,j,n,dtsur,dxsur,tmp},
If[!(Dimensions[sur1]==
Dimensions[sur2]==
Dimensions[sur3]),
Message[ParametricListPlot3D::wrdim];Return[$Failed];];
{nxpoints,ntpoints}=Dimensions[sur1];
sur[1]=sur1;sur[2]=sur2;sur[3]=sur3;nopnop;
Do[ dtsur[i]=Table[ND[sur[i][[n]],tdelta],{n,nxpoints}],
{i,3}];
Do[
sur[i]=Transpose[sur[i]];
dxsur[i]=Table[ND[sur[i][[j]],xdelta],{j,ntpoints}];
dxsur[i]=Transpose[dxsur[i]],
{i,3}];
gg2=Table[ng[n]^2,{n,nxpoints}];
nop;
{ListPlot3D[Evaluate[
(ep*dtsur[1]^2+dtsur[2]^2+dtsur[3]^2)-ep
],PlotLabel->"Errore su E.\n\n\n\n\n",opts],
ListPlot3D[Evaluate[
(ep*dtsur[1] * dxsur[1] +
dtsur[2] * dxsur[1] +
dtsur[3] * dxsur[1])
],PlotLabel->"Errore su F.\n\n\n\n\n",opts],
ListPlot3D[Evaluate[
(ep*dxsur[1]^2+dxsur[2]^2+dxsur[3]^2)-gg2
],PlotLabel->"Errore su G.\n\n\n\n\n",opts]}
];
(*
:[font = input; initialization; preserveAspect]
*)
ImmergiEdisegna[ep_?EpQ,g_,
{t_Symbol,tmin_?NumberQ,tmax_?NumberQ},ntpoints_Integer,
{x_Symbol,xmin_?NumberQ,xmax_?NumberQ},nxpoints_Integer,
gradoint_Integer,nome_String]:=
Block[{stretch,rc,k,opts,cck,numk,sur,ng,err,gopts,
sup,leg,col0,th0,th1,th2,y,tdelta,xdelta,maxk,mink},
Print["-> Calcolo gli angoli theta."];
opts=Sequence[SamplePoints$t->ntpoints,
SamplePoints$x->nxpoints
];
rc=ImmersionOf[g,{t,tmin,tmax},{x,xmin,xmax},
{th0,th1,th2},ng,{tdelta,xdelta},ep,opts];
If[rc==$Failed,Print["Errore nel calcolo degli angoli"];Return[]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Calcolo la superfice "];
rc=ImmersedSurface[{th0,th1,th2},ng,
{tdelta,xdelta},ep,gradoint,sur,opts];
If[rc==$Failed,Print["Errore nel calcolo eq. Superfice"];Return[]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Calcolo la curvatura."];
k=Kurvature[1,g,{t,x},ep];
cck=Compile[{t,x},k];
numk=Table[cck[t,x],{x,xmin,xmax,xdelta},
{t,tmin,tmax,tdelta}];
maxk=Max[numk];mink=Min[numk];
If[maxk==0. && mink==0., maxk=0.1;mink=-0.1];
col0=(maxk+mink)/2;
stretch=(maxk-col0);
If[Chop[maxk-mink]==0,stretch=1;maxk=maxk*1.1;mink=mink*0.9,
stretch=1.5/stretch;];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Disegno e coloro la superfice"];
sup=ParametricListPlot3D[sur[1],sur[2],sur[3],
Map[(ColorStretched[#-col0,stretch])&,
numk,{2}],PlotLabel->"E="<>ToString[ep]<>
" G(t,x)="<>
ToString[InputForm[g^2]]<>"\n\n\n",
AxesLabel->{ "\nx0",
"x1 \n\n",
"x2 "},
AxesEdge->{ {-1,-1},
{-1,1},
{-1,-1}},
Axes->True];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Salvo la superfice -sup- in ",Evaluate[nome]];
Put[Definition["sup"], Evaluate[nome]];
Print["-> Disegno la legenda -leg- dei colori e salvo."];
leg=PlotLegend[mink,maxk,col0,stretch];
PutAppend[Definition["leg"], Evaluate[nome]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Disegno l'errore -err- commesso e salvo."];
err=PlotImmersionError[sur[1],sur[2],sur[3],
ng,ep,{tdelta,xdelta},
AxesLabel->{ "j ",
"n ",
"Err "},
AxesEdge->{ {-1,-1},
{-1,1},
{-1,-1}}];
PutAppend[Definition["err"],Evaluate[nome]];
PutAppend[Definition["g","k","tdelta",
"xdelta","tmin","xmin",
"tmax","xmax"], Evaluate[nome]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["********FINE***********"];
];
(*
;[s]
1:0,0;2782,-1;
1:1,0,0 ,Courier,0,10,0,0,0;
:[font = input; initialization; preserveAspect]
*)
ImmergiEdisegnaBatch[ep_?EpQ,g_,
{t_Symbol,tmin_?NumberQ,tmax_?NumberQ},ntpoints_Integer,
{x_Symbol,xmin_?NumberQ,xmax_?NumberQ},nxpoints_Integer,
gradoint_Integer,nome_String]:=
Block[{stretch,rc,k,opts,cck,numk,sur,ng,err,gopts,
sup,leg,col0,th0,th1,th2,y,tdelta,xdelta,maxk,mink},
Print["-> Calcolo gli angoli theta."];
opts=Sequence[SamplePoints$t->ntpoints,
SamplePoints$x->nxpoints
];
rc=ImmersionOf[g,{t,tmin,tmax},{x,xmin,xmax},
{th0,th1,th2},ng,{tdelta,xdelta},ep,opts];
If[rc==$Failed,Print["Errore nel calcolo degli angoli"];Return[]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Calcolo la superfice "];
rc=ImmersedSurface[{th0,th1,th2},ng,
{tdelta,xdelta},ep,gradoint,sur,opts];
If[rc==$Failed,Print["Errore nel calcolo eq. Superfice"];Return[]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Calcolo la curvatura."];
k=Kurvature[1,g,{t,x},ep];
cck=Compile[{t,x},k];
numk=Table[cck[t,x],{x,xmin,xmax,xdelta},
{t,tmin,tmax,tdelta}];
maxk=Max[numk];mink=Min[numk];
If[maxk==0. && mink==0., maxk=0.1;mink=-0.1];
col0=(maxk+mink)/2;
stretch=(maxk-col0);
If[Chop[maxk-mink]==0,stretch=1;maxk=maxk*1.1;mink=mink*0.9,
stretch=1.5/stretch;];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Disegno e coloro la superfice"];
sup=ParametricListPlot3D[sur[1],sur[2],sur[3],
Map[(ColorStretched[#-col0,stretch])&,
numk,{2}],
AxesLabel->{ "\nx0",
"x1 \n\n",
"x2 "},
AxesEdge->{ {-1,-1},
{-1,1},
{-1,-1}},
Axes->True,
DefaultFont->{"Symbol",10},
DisplayFunction->Function[
Display[nome<>"_sup.ps", #]]
];
RPLScript[nome<>"_sup.rpl",sup];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Salvo la superfice -sup- in ",Evaluate[nome]];
Put[Definition["sup"], Evaluate[nome]];
Print["-> Disegno la legenda -leg- dei colori e salvo."];
leg=PlotLegend[mink,maxk,col0,stretch];
Show[leg,DisplayFunction->Function[Display[nome<>"_leg.ps", #]]];
PutAppend[Definition["leg"], Evaluate[nome]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["-> Disegno l'errore -err- commesso e salvo."];
err=PlotImmersionError[sur[1],sur[2],sur[3],
ng,ep,{tdelta,xdelta},
AxesLabel->{ "j ",
"n ",
"Err "},
AxesEdge->{ {-1,-1},
{-1,1},
{-1,-1}},
DisplayFunction->Identity];
Show[err[[1]],DisplayFunction->Function[Display[nome<>"_eer.ps", #]]];
Show[err[[2]],DisplayFunction->Function[Display[nome<>"_fer.ps", #]]];
Show[err[[3]],DisplayFunction->Function[Display[nome<>"_ger.ps", #]]];
PutAppend[Definition["err"],Evaluate[nome]];
PutAppend[Definition["g","k","tdelta",
"xdelta","tmin","xmin",
"tmax","xmax"], Evaluate[nome]];
Print["~> Ottimizzo la memoria."];
Print[" Guadagnati... ",nop," bytes."];
Print["********FINE***********"];
];
(*
;[s]
3:0,0;1816,1;1817,2;3139,-1;
3:1,0,0 ,Courier,0,10,0,0,0;1,0,0 ,Courier,1,12,0,0,0;1,0,0 ,Courier,0,10,0,0,0;
^*)