/* This file is part of the source code for 3D-XplorMath-J, Version 1.0 (January 2008).
* Copyright (c) 2008 The 3D-XplorMath Consortium (http://3d-xplormath.org).
* This source code is released under a BSD License, which allows redistribution
* in source and binary form, with or without modification, provided copyright
* and license information are included, and with no warranty or guarantee of
* any kind. For details, see http://3d-xplormath.org/j/source/BSDLicense.txt
*/
package vmm3d.spacecurve.parametric;
import vmm3d.actions.ActionList;
import vmm3d.actions.ActionRadioGroup;
import vmm3d.core.I18n;
import vmm3d.core.Parameter;
import vmm3d.core.RealParamAnimateable;
import vmm3d.core.View;
import vmm3d.core3D.Vector3D;
/**
* Defines a space curve of constant curvature kappa by integrating the Frenet equations with kappa=aa
* and a torsion function tau(t) = bb + cc*sin(t) + dd*sin(2t) + ee*sin(3t).
* This Class behaves like an explicitly parametrized curve, because the program precomputes the
* Frenet frame on 100 points between tmin and tmax so that, when the Superclass calls the method
* Vector3D value(), the ODEsolver only has to integrate from the nearest precomputed point.
*/
public class ConstantCurvature extends SpaceCurveParametric {
RealParamAnimateable aa = new RealParamAnimateable("vmm.spacecurve.parametric.ConstantCurvature.aa",1.0, 1.0, 1.0);
RealParamAnimateable bb = new RealParamAnimateable("vmm.spacecurve.parametric.ConstantCurvature.bb",0.0, 0.0, 0.0);
RealParamAnimateable cc = new RealParamAnimateable("vmm.spacecurve.parametric.ConstantCurvature.cc",1.00844728,1.354743, 0.314);
RealParamAnimateable dd = new RealParamAnimateable("vmm.spacecurve.parametric.ConstantCurvature.dd",0.0, 0.0, 0.0);
RealParamAnimateable ee = new RealParamAnimateable("vmm.spacecurve.parametric.ConstantCurvature.ee",0.0, 0.0, 0.0);
double kappa = 0.0;
double b = 0.0;
double c = 0.0;
double d = 0.0;
double e = 0.0;
private int exampleNumber = 0;
private ActionRadioGroup exampleSelect;
public ConstantCurvature() {
setDefaultViewpoint(new Vector3D(0,-16,-6));
setDefaultWindow(-4,4,-3,3);
tResolution.setValueAndDefault(300);
tmin.setValueAndDefault(0);
tmax.setValueAndDefaultFromString("4*pi");
addParameter(ee);
addParameter(dd);
addParameter(cc);
addParameter(bb);
addParameter(aa);
tubeSize.setValueAndDefault(0.2);
// Next the constructor part for the radio buttons:
ActionRadioGroup group = new ActionRadioGroup() {
public void optionSelected(int selectedIndex) {
setExampleNumberFunction(selectedIndex);
}
};
group.addItem(I18n.tr("vmm.spacecurve.parametric.ConstantCurvature.example0"));
group.addItem(I18n.tr("vmm.spacecurve.parametric.ConstantCurvature.example1"));
group.addItem(I18n.tr("vmm.spacecurve.parametric.ConstantCurvature.example2"));
group.addItem(I18n.tr("vmm.spacecurve.parametric.ConstantCurvature.example3"));
group.addItem(I18n.tr("vmm.spacecurve.parametric.ConstantCurvature.example4"));
group.setSelectedIndex(0);
exampleSelect = group;
}
// First the code connected with the radio buttons:
public int getExampleNumberFunction() {
return exampleNumber;
}
public void setExampleNumberFunction(int exampleNumber) {
if (this.exampleNumber != exampleNumber) {
this.exampleNumber = exampleNumber;
//exampleSelect.setSelectedIndex(exampleNumber);
System.out.println("exampleNumber = "+exampleNumber);
resetFourierConstants(exampleNumber); // here we need to prepare the selected example
forceRedraw();
}
}
public ActionList getActionsForView(final View view) {
ActionList actions = super.getActionsForView(view);
actions.add(null); // separator in menu
ActionList submenu = new ActionList(I18n.tr("vmm.spacecurve.parametric.ConstantCurvature.closedCurves"));
submenu.add(exampleSelect);
actions.add(submenu);
return actions;
}
private void resetFourierConstants(int examplNum){
needsNewArray = true;
tResolution.reset(300);
aa.reset(1.0);
switch (examplNum){
case 0:
tmax.reset("4*pi");
bb.reset(0.0);
cc.reset(1.0084472);
break;
case 1:
tmax.reset("6*pi");
bb.reset(0.0);
cc.reset(1.3547);
break;
case 2:
tmax.reset("10*pi");
bb.reset(0.2470);
cc.reset(0.4053);
break;
case 3:
tmax.reset("10*pi");
bb.reset(0.5117);
cc.reset(0.8551);
break;
case 4:
tResolution.reset(600);
tmax.reset("22*pi");
aa.reset(0.8944);
bb.reset(0.5005);
cc.reset(0.4129);
break;
}
getConstants();
}
// The code for the curve starts here
private void getConstants(){
tr = tResolution.getValue();
kappa = aa.getValue();
b = bb.getValue();
c = cc.getValue();
d = dd.getValue();
e = ee.getValue();
}
int tr = tResolution.getValue();
Vector3D startPoint = new Vector3D(0,0,0);
Vector3D P = new Vector3D(0,0,-1);
Vector3D T = new Vector3D(1,0,0);
Vector3D N = new Vector3D(0,1,0);
Vector3D B = new Vector3D(0,0,1);
Vector3D[] initialFrenet = new Vector3D [] {startPoint,T,N,B};
Vector3D[] repere = new Vector3D[]{P, T, N, B};
Vector3D[] pointSet = new Vector3D[tr+1];
Vector3D[] helperFrame = new Vector3D[4*(tr + 1)];
boolean needsNewArray = true;
/*{
System.out.println(tr);
}*/
public void parameterChanged(Parameter param, Object oldValue, Object newValue) {
super.parameterChanged(param, oldValue, newValue);
needsNewArray = true;
}
protected double tau(double t){
return (b + c*Math.sin(t) + d*Math.sin(2*t) + e*Math.sin(3*t));
}
/*protected double dtau(double t){
return (c*Math.cos(t) + 2*d*Math.cos(2*t) + 3*e*Math.cos(3*t));
}*/
/**
* Second order method, "Step to midpoint"
*/
protected Vector3D[] frenetODEstep2(double tInitial, double tFinal, Vector3D[]initialVal, int numSubdivision){
P = initialVal[0];
T = initialVal[1];
N = initialVal[2];
B = initialVal[3];
double t = tInitial;
double newSub = (1.0+Math.abs(tFinal-tInitial))*numSubdivision;
double dt = (tFinal-tInitial)/newSub;
Vector3D midT = new Vector3D(T);
Vector3D midN = new Vector3D(N);
Vector3D midB = new Vector3D(B);
for (int i = 0; i < newSub ; i++)
{
midT = T.plus(N.times(kappa*dt/2));
midN = N.plus(T.times(-kappa*dt/2)).plus(B.times(tau(t)*dt/2));
midB = B.plus(N.times(-tau(t)*dt/2));
P = P.plus( T.plus(midT.times(4)).times(dt/6.0) );
t = t+dt/2;
T = T.plus(midN.times(kappa*dt));
N = N.plus(midT.times(-kappa*dt)).plus(midB.times(tau(t)*dt));
B = B.plus(midN.times(-tau(t)*dt));
P = P.plus(T.times(dt/6)); // Simpson
t=t+dt/2;
}
return new Vector3D[]{P,T,N,B};
}
/**
* Runge Kutta integrator, 4th order
*/
Vector3D dT1, dN1, dB1, dT2, dN2, dB2, dT3, dN3, dB3, dT4, dN4, dB4, Vvv;
{dT1=new Vector3D(0,0,0); dN1=new Vector3D(0,0,0); dB1=new Vector3D(0,0,0);
dT2=new Vector3D(0,0,0); dN2=new Vector3D(0,0,0); dB2=new Vector3D(0,0,0);
dT3=new Vector3D(0,0,0); dN3=new Vector3D(0,0,0); dB3=new Vector3D(0,0,0);
dT4=new Vector3D(0,0,0); dN4=new Vector3D(0,0,0); dB4=new Vector3D(0,0,0);
Vvv=new Vector3D(0,0,0);}
protected Vector3D[] frenetODEstep4(double tInitial, double tFinal, Vector3D[]initialVal, int numSubdivision){
P = initialVal[0];
T = initialVal[1];
N = initialVal[2];
B = initialVal[3];
double t = tInitial;
double newSub = Math.floor( (1.0+Math.abs(tFinal-tInitial))*numSubdivision );
double dt = (tFinal-tInitial)/newSub/2.0;
newSub = newSub - 1.0/16.0;
for (int i = 0; i < newSub; i++) {
for (int j = 1; j < 3; j++) { // Two Runge-Kutta steps for one Simpson 1-4-1
dT1.assignTimes(kappa*dt/2, N);
dN1.assignLinComb(-kappa*dt/2, T, tau(t)*dt/2, B); //T.times(-kappa*dt/2).plus(B.times(tau(t)*dt/2));
dB1.assignTimes(-tau(t)*dt/2, N); //N.times(-tau(t)*dt/2);
t = t+dt/2;
dT2.assignSumTimes(N, dN1, kappa*dt); //(N.plus(dN1)).times(kappa*dt);
dN2.assignSumTimes(T, dT1, -kappa*dt); Vvv.assignSumTimes(B, dB1, tau(t)*dt);
dN2.assignPlus(Vvv); //(T.plus(dT1)).times(-kappa*dt).plus((B.plus(dB1)).times(tau(t)*dt));
dB2.assignSumTimes(N, dN1, -tau(t)*dt); //(N.plus(dN1)).times(-tau(t)*dt);
dT3.assignLinComb(kappa*dt,N, kappa*dt/2,dN2); //(N.plus(dN2.times(0.5))).times(kappa*dt);
dN3.assignLinComb(-kappa*dt,T, -kappa*dt/2,dT2); Vvv.assignLinComb(tau(t)*dt,B, tau(t)*dt/2,dB2);
dN3.assignPlus(Vvv); //(T.plus(dT2.times(0.5))).times(-kappa*dt).plus((B.plus(dB2.times(0.5))).times(tau(t)*dt));
dB3.assignLinComb(-tau(t)*dt,N, -tau(t)*dt/2,dN2); //(N.plus(dN2.times(0.5))).times(-tau(t)*dt);
t=t+dt/2;
dT4.assignSumTimes(N, dN3, kappa*dt/2); //(N.plus(dN3)).times(kappa*dt/2);
dN4.assignSumTimes(T, dT3, -kappa*dt/2); Vvv.assignSumTimes(B, dB3, tau(t)*dt/2);
dN4.assignPlus(Vvv); //(T.plus(dT3)).times(-kappa*dt/2).plus((B.plus(dB3)).times(tau(t)*dt/2));
dB4.assignSumTimes(N, dN3, -tau(t)*dt/2); //(N.plus(dN3)).times(-tau(t)*dt/2);
// now one needs to return new vectors:
P = P.plus(T.times(dt*j*j/3.0));// part of Simpson .assignLinComb(1, j*j*dt/3.0, T); //
Vvv.assignSumTimes(dT1, dT2, dT3, dT4, 1.0/3.0);
T = T.plus(Vvv);
// T = T.plus( ( dT1.plus(dT2.plus(dT3.plus(dT4))) ).times(1.0/3.0) );
Vvv.assignSumTimes(dN1, dN2, dN3, dN4, 1.0/3.0);
N = N.plus(Vvv);
// N = N.plus( ( dN1.plus(dN2.plus(dN3.plus(dN4))) ).times(1.0/3.0) );
Vvv.assignSumTimes(dB1, dB2, dB3, dB4, 1.0/3.0);
B = B.plus(Vvv);
// B = B.plus( ( dB1.plus(dB2.plus(dB3.plus(dB4))) ).times(1.0/3.0) );
}
P.assignLinComb(1, dt/3.0, T);
//P.plus(T.times(dt/3.0)); // Simpson integration. Because of the two Runge Kutta steps, dt is half the step size
}
return new Vector3D[]{P,T,N,B};
}
private void createFrenetArray(){
getConstants();
double dt = (tmax.getValue() - tmin.getValue())/tr;
double t = tmin.getValue();
helperFrame = new Vector3D [4*(tr + 1)];
pointSet = new Vector3D[tr+1];
helperFrame[0]= new Vector3D(0,0,0);
helperFrame[1]= new Vector3D(1,0,0);
helperFrame[2]= new Vector3D(0,1,0);
helperFrame[3]= new Vector3D(0,0,1);
pointSet[0] = new Vector3D(helperFrame[0]);
for (int i = 0; i < tr ; i++) { //tr = tResolution.getValue()
for (int j = 0; j < 4 ; j++){
initialFrenet[j] = helperFrame[j+4*i];
}
repere = frenetODEstep4(t,t+dt,initialFrenet,32);
for (int j = 0; j < 4 ; j++){
helperFrame[j+4*(i+1)] = new Vector3D(repere[j]);
}
pointSet[i+1] = (helperFrame[4*(i+1)]);
t = t+dt;
}
startPoint = getCenterOfPoints(pointSet, tr+1);
startPoint.x = - startPoint.x; startPoint.y = - startPoint.y; startPoint.z = - startPoint.z;
for (int i = 0; i <= tr ; i++) {
helperFrame[4*i].assignPlus(startPoint);
}
// See the line 'needsNewArray = true;' in the above override of public void parameterChanged(...
}
// The following override method does not use the default numerical differentiation but integrates
// the Frenet ODE. The precomputed array is used to speed up the computation.
public Vector3D[] makeRepereMobile(double t) {
int i = 0;
double tStart = 0.0;
if (needsNewArray) {
createFrenetArray();
needsNewArray = false;
}
i = (int) Math.floor(tResolution.getValue()*(t-tmin.getValue())/(tmax.getValue()-tmin.getValue()));
i = Math.max(0, Math.min(i, tResolution.getValue()));
tStart = tmin.getValue() + i*(tmax.getValue()-tmin.getValue())/tResolution.getValue();
for (int j = 0; j < 4 ; j++){
initialFrenet[j] = helperFrame[j+ 4*i];
}
//repere = frenetODEstep2(tStart,t,initialFrenet,128);
repere = frenetODEstep4(tStart,t,initialFrenet,16);
return repere;
}
// This method has to be supplied by the Class. It behaves like a function call,
// but the values are produced by solving an ODE.
protected Vector3D value(double t) {
int i = 0;
double tStart = 0.0;
if (needsNewArray) {
createFrenetArray();
needsNewArray = false;
}
i = (int) Math.floor(tResolution.getValue()*(t-tmin.getValue())/(tmax.getValue()-tmin.getValue()));
i = Math.max(0, Math.min(i, tResolution.getValue()));
tStart = tmin.getValue() + i*(tmax.getValue()-tmin.getValue())/tResolution.getValue();
for (int j = 0; j < 4 ; j++){
initialFrenet[j] = helperFrame[j+ 4*i];
}
//Vector3D point=new Vector3D(frenetODEstep2(tStart,t,initialFrenet,128)[0]);
Vector3D point=(frenetODEstep4(tStart,t,initialFrenet,16)[0]);
return point;
}
}
Ads help maintain this website.