Matlab-Programs/SaveSatePythonTest.py

+# state file generated using paraview version 5.7.0
+
+# ----------------------------------------------------------------
+# setup views used in the visualization
+# ----------------------------------------------------------------
+
+# trace generated using paraview version 5.7.0
+#
+# To ensure correct image size when batch processing, please search 
+# for and uncomment the line `# renderView*.ViewSize = [*,*]`
+
+#### import the simple module from the paraview
+from paraview.simple import *
+#### disable automatic camera reset on 'Show'
+paraview.simple._DisableFirstRenderCameraReset()
+
+# Create a new 'Render View'
+renderView1 = CreateView('RenderView')
+renderView1.ViewSize = [964, 795]
+renderView1.AxesGrid = 'GridAxes3DActor'
+renderView1.CenterOfRotation = [5.5, 5.5, 5.5]
+renderView1.StereoType = 'Crystal Eyes'
+renderView1.CameraPosition = [-1.9013843868480158, -13.320979633549339, 35.98066109404487]
+renderView1.CameraFocalPoint = [1.538094468601466, -4.574729188366619, 21.81606919426716]
+renderView1.CameraViewUp = [0.05745797827794325, -0.8556355497369458, -0.5143796134748793]
+renderView1.CameraFocalDisk = 1.0
+renderView1.CameraParallelScale = 7.794228634059948
+renderView1.Background = [0.32, 0.34, 0.43]
+
+SetActiveView(None)
+
+# ----------------------------------------------------------------
+# setup view layouts
+# ----------------------------------------------------------------
+
+# create new layout object 'Layout #1'
+layout1 = CreateLayout(name='Layout #1')
+layout1.AssignView(0, renderView1)
+
+# ----------------------------------------------------------------
+# restore active view
+SetActiveView(renderView1)
+# ----------------------------------------------------------------
+
+# ----------------------------------------------------------------
+# setup the data processing pipelines
+# ----------------------------------------------------------------
+
+# create a new 'Legacy VTK Reader'
+windvtk = LegacyVTKReader(FileNames=['/home/klaus/Desktop/DUNE/dune-microstructure/Matlab-Programs/wind.vtk'])
+
+# ----------------------------------------------------------------
+# setup the visualization in view 'renderView1'
+# ----------------------------------------------------------------
+
+# show data from windvtk
+windvtkDisplay = Show(windvtk, renderView1)
+
+# get color transfer function/color map for 'type'
+typeLUT = GetColorTransferFunction('type')
+typeLUT.RGBPoints = [1.0, 0.0, 0.0, 1.0, 1.0, 0.17254901960784313, 0.043137254901960784, 0.9215686274509803, 1.0, 0.785929, 0.247056, 0.295477, 1.0, 0.17254901960784313, 0.043137254901960784, 0.9215686274509803, 1.1698113679885864, 0.41568627450980394, 0.00392156862745098, 0.9686274509803922, 1.2672955989837646, 0.48627450980392156, 0.0, 0.9725490196078431, 1.4937106370925903, 1.0, 0.058823529411764705, 0.984313725490196, 1.7672955989837646, 1.0, 0.043137254901960784, 0.5529411764705883, 2.0, 1.0, 0.0, 0.01568627450980392]
+typeLUT.NanColor = [0.0, 1.0, 0.0]
+typeLUT.NumberOfTableValues = 593
+typeLUT.ScalarRangeInitialized = 1.0
+
+# get opacity transfer function/opacity map for 'type'
+typePWF = GetOpacityTransferFunction('type')
+typePWF.Points = [1.0, 1.0, 0.5, 0.0, 1.2830188274383545, 1.0, 0.5, 0.0, 1.4654088020324707, 1.0, 0.5, 0.0, 1.6289308071136475, 1.0, 0.5, 0.0, 1.7987420558929443, 0.9926470518112183, 0.5, 0.0, 2.0, 1.0, 0.5, 0.0]
+typePWF.ScalarRangeInitialized = 1
+
+# trace defaults for the display properties.
+windvtkDisplay.Representation = 'Points'
+windvtkDisplay.ColorArrayName = ['POINTS', 'type']
+windvtkDisplay.LookupTable = typeLUT
+windvtkDisplay.PointSize = 3.0
+windvtkDisplay.OSPRayScaleArray = 'type'
+windvtkDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
+windvtkDisplay.SelectOrientationVectors = 'None'
+windvtkDisplay.ScaleFactor = 0.9
+windvtkDisplay.SelectScaleArray = 'type'
+windvtkDisplay.GlyphType = 'Arrow'
+windvtkDisplay.GlyphTableIndexArray = 'type'
+windvtkDisplay.GaussianRadius = 0.045
+windvtkDisplay.SetScaleArray = ['POINTS', 'type']
+windvtkDisplay.ScaleTransferFunction = 'PiecewiseFunction'
+windvtkDisplay.OpacityArray = ['POINTS', 'type']
+windvtkDisplay.OpacityTransferFunction = 'PiecewiseFunction'
+windvtkDisplay.DataAxesGrid = 'GridAxesRepresentation'
+windvtkDisplay.PolarAxes = 'PolarAxesRepresentation'
+windvtkDisplay.ScalarOpacityFunction = typePWF
+windvtkDisplay.ScalarOpacityUnitDistance = 0.08660254037844385
+
+# init the 'PiecewiseFunction' selected for 'OSPRayScaleFunction'
+windvtkDisplay.OSPRayScaleFunction.Points = [1.0, 0.8602941036224365, 0.5, 0.0, 1.213836476688135, 0.625, 0.5, 0.0, 1.3616352081298828, 0.6691176295280457, 0.5, 0.0, 1.6666333299996667, 0.8676470518112183, 0.5, 0.0, 1.7358490228652954, 0.6911764740943909, 0.5, 0.0, 2.0, 0.8014705777168274, 0.5, 0.0]
+
+# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
+windvtkDisplay.ScaleTransferFunction.Points = [1.0, 0.8602941036224365, 0.5, 0.0, 1.213836476688135, 0.625, 0.5, 0.0, 1.3616352081298828, 0.6691176295280457, 0.5, 0.0, 1.6666333299996667, 0.8676470518112183, 0.5, 0.0, 1.7358490228652954, 0.6911764740943909, 0.5, 0.0, 2.0, 0.8014705777168274, 0.5, 0.0]
+
+# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
+windvtkDisplay.OpacityTransferFunction.Points = [1.0, 0.8602941036224365, 0.5, 0.0, 1.213836476688135, 0.625, 0.5, 0.0, 1.3616352081298828, 0.6691176295280457, 0.5, 0.0, 1.6666333299996667, 0.8676470518112183, 0.5, 0.0, 1.7358490228652954, 0.6911764740943909, 0.5, 0.0, 2.0, 0.8014705777168274, 0.5, 0.0]
+
+# setup the color legend parameters for each legend in this view
+
+# get color legend/bar for typeLUT in view renderView1
+typeLUTColorBar = GetScalarBar(typeLUT, renderView1)
+typeLUTColorBar.Position = [0.8796680497925311, 0.01509433962264151]
+typeLUTColorBar.Title = 'type'
+typeLUTColorBar.ComponentTitle = ''
+
+# set color bar visibility
+typeLUTColorBar.Visibility = 1
+
+# show color legend
+windvtkDisplay.SetScalarBarVisibility(renderView1, True)
+
+# ----------------------------------------------------------------
+# setup color maps and opacity mapes used in the visualization
+# note: the Get..() functions create a new object, if needed
+# ----------------------------------------------------------------
+
+# ----------------------------------------------------------------
+# finally, restore active source
+SetActiveSource(windvtk)
+# ----------------------------------------------------------------
\ No newline at end of file

Matlab-Programs/Stefan_parabolic.m

+
+clc
+q1=1;
+q2=2;
+q12=1/2;
+q3=((4*q1*q2)^(1/2)-q12)/2;
+H=[2*q1,q12+2*q3;q12+2*q3,2*q2];
+A=[q1,1/2*q12;1/2*q12,q2];
+abar=[q12+2*q3;2*q2];
+abar=(abar(1)^2+abar(2)^2)^(-1/2).*abar
+b=1*A\abar
+sstar=1/(q1+q2)*abar'*(A*b)
+abarperp=[abar(2);-abar(1)]
+%
+N=10;
+T=linspace(-sstar*(q12+2*q3)/(2*q2),sstar*(2*q2)/(q12+2*q3),N)
+abars=sstar.*abar+T.*abarperp
+%
+G=[];
+e=[];
+alpha=[];
+kappa=[];
+K=[];
+for k=1:N
+    if abars(1,k)*abars(2,k)>=0
+        K=[K,k];
+        G=[G,[abars(1,k);abars(2,k);sqrt(2*abars(1,k)*abars(2,k))]];
+        g=[sqrt(abars(1,k)/(abars(1,k)+abars(2,k)));sqrt(abars(2,k)/(abars(1,k)+abars(2,k)))];
+        e=[e,g];
+        kappa=[kappa,abars(1,k)+abars(2,k)];
+%         alpha=[alpha,atan2(g(1),g(2))];                                                 %reversed here ???? 
+        alpha=[alpha,atan2(g(2),g(1))];    
+    end
+end
+% G1=[1,0;0,0];G2=[0,0;0,1];G3=[0,1/sqrt(2);1/sqrt(2),0];
+alpha
+kappa
+min(alpha)
+max(alpha)
+hold off;
+plot(alpha,kappa);
+hold;
+plot(-alpha,kappa);

Matlab-Programs/T.py

+# state file generated using paraview version 5.7.0
+
+# ----------------------------------------------------------------
+# setup views used in the visualization
+# ----------------------------------------------------------------
+
+# trace generated using paraview version 5.7.0
+#
+# To ensure correct image size when batch processing, please search 
+# for and uncomment the line `# renderView*.ViewSize = [*,*]`
+
+#### import the simple module from the paraview
+from paraview.simple import *
+#### disable automatic camera reset on 'Show'
+paraview.simple._DisableFirstRenderCameraReset()
+
+# Create a new 'Render View'
+renderView1 = CreateView('RenderView')
+renderView1.ViewSize = [901, 795]
+renderView1.AxesGrid = 'GridAxes3DActor'
+renderView1.CenterOfRotation = [0.22986745834350586, 0.03941011428833008, 0.17028677463531494]
+renderView1.StereoType = 'Crystal Eyes'
+renderView1.CameraPosition = [53.44770716226582, -21.61022632014142, -4.102774671003825]
+renderView1.CameraFocalPoint = [0.22986745834350583, 0.039410114288330106, 0.17028677463531497]
+renderView1.CameraViewUp = [-0.28235039193869665, -0.7988816556360265, 0.531099196441028]
+renderView1.CameraFocalDisk = 1.0
+renderView1.CameraParallelScale = 5.748834779828519
+renderView1.Background = [0.32, 0.34, 0.43]
+
+SetActiveView(None)
+
+# ----------------------------------------------------------------
+# setup view layouts
+# ----------------------------------------------------------------
+
+# create new layout object 'Layout #1'
+layout1 = CreateLayout(name='Layout #1')
+layout1.AssignView(0, renderView1)
+
+# ----------------------------------------------------------------
+# restore active view
+SetActiveView(renderView1)
+# ----------------------------------------------------------------
+
+# ----------------------------------------------------------------
+# setup the data processing pipelines
+# ----------------------------------------------------------------
+
+# create a new 'Legacy VTK Reader'
+matlab_exportvtk = LegacyVTKReader(FileNames=['/home/klaus/Desktop/DUNE/dune-microstructure/Matlab-Programs/matlab_export.vtk'])
+
+# create a new 'Legacy VTK Reader'
+pointSetvtk = LegacyVTKReader(FileNames=['/home/klaus/Desktop/Paraview-Data/ParaView-v5.9.0/Testing/Data/pointSet.vtk'])
+
+# ----------------------------------------------------------------
+# setup the visualization in view 'renderView1'
+# ----------------------------------------------------------------
+
+# show data from pointSetvtk
+pointSetvtkDisplay = Show(pointSetvtk, renderView1)
+
+# get color transfer function/color map for 'RandomPointScalars'
+randomPointScalarsLUT = GetColorTransferFunction('RandomPointScalars')
+randomPointScalarsLUT.RGBPoints = [2.0, 0.231373, 0.298039, 0.752941, 5.5, 0.865003, 0.865003, 0.865003, 9.0, 0.705882, 0.0156863, 0.14902]
+randomPointScalarsLUT.ScalarRangeInitialized = 1.0
+
+# trace defaults for the display properties.
+pointSetvtkDisplay.Representation = 'Points'
+pointSetvtkDisplay.ColorArrayName = ['POINTS', 'RandomPointScalars']
+pointSetvtkDisplay.LookupTable = randomPointScalarsLUT
+pointSetvtkDisplay.RenderPointsAsSpheres = 1
+pointSetvtkDisplay.OSPRayScaleArray = 'RandomPointScalars'
+pointSetvtkDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
+pointSetvtkDisplay.SelectOrientationVectors = 'RandomCellVectors'
+pointSetvtkDisplay.ScaleFactor = 0.7009325742721558
+pointSetvtkDisplay.SelectScaleArray = 'RandomPointScalars'
+pointSetvtkDisplay.GlyphType = 'Arrow'
+pointSetvtkDisplay.GlyphTableIndexArray = 'RandomPointScalars'
+pointSetvtkDisplay.GaussianRadius = 0.03504662871360779
+pointSetvtkDisplay.SetScaleArray = ['POINTS', 'RandomPointScalars']
+pointSetvtkDisplay.ScaleTransferFunction = 'PiecewiseFunction'
+pointSetvtkDisplay.OpacityArray = ['POINTS', 'RandomPointScalars']
+pointSetvtkDisplay.OpacityTransferFunction = 'PiecewiseFunction'
+pointSetvtkDisplay.DataAxesGrid = 'GridAxesRepresentation'
+pointSetvtkDisplay.PolarAxes = 'PolarAxesRepresentation'
+
+# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
+pointSetvtkDisplay.ScaleTransferFunction.Points = [2.0, 0.0, 0.5, 0.0, 9.0, 1.0, 0.5, 0.0]
+
+# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
+pointSetvtkDisplay.OpacityTransferFunction.Points = [2.0, 0.0, 0.5, 0.0, 9.0, 1.0, 0.5, 0.0]
+
+# show data from matlab_exportvtk
+matlab_exportvtkDisplay = Show(matlab_exportvtk, renderView1)
+
+# trace defaults for the display properties.
+matlab_exportvtkDisplay.Representation = 'Surface With Edges'
+matlab_exportvtkDisplay.ColorArrayName = [None, '']
+matlab_exportvtkDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
+matlab_exportvtkDisplay.SelectOrientationVectors = 'None'
+matlab_exportvtkDisplay.ScaleFactor = 9.9
+matlab_exportvtkDisplay.SelectScaleArray = 'None'
+matlab_exportvtkDisplay.GlyphType = 'Arrow'
+matlab_exportvtkDisplay.GlyphTableIndexArray = 'None'
+matlab_exportvtkDisplay.GaussianRadius = 0.495
+matlab_exportvtkDisplay.SetScaleArray = [None, '']
+matlab_exportvtkDisplay.ScaleTransferFunction = 'PiecewiseFunction'
+matlab_exportvtkDisplay.OpacityArray = [None, '']
+matlab_exportvtkDisplay.OpacityTransferFunction = 'PiecewiseFunction'
+matlab_exportvtkDisplay.DataAxesGrid = 'GridAxesRepresentation'
+matlab_exportvtkDisplay.PolarAxes = 'PolarAxesRepresentation'
+
+# setup the color legend parameters for each legend in this view
+
+# get color legend/bar for randomPointScalarsLUT in view renderView1
+randomPointScalarsLUTColorBar = GetScalarBar(randomPointScalarsLUT, renderView1)
+randomPointScalarsLUTColorBar.Title = 'RandomPointScalars'
+randomPointScalarsLUTColorBar.ComponentTitle = ''
+
+# set color bar visibility
+randomPointScalarsLUTColorBar.Visibility = 1
+
+# show color legend
+pointSetvtkDisplay.SetScalarBarVisibility(renderView1, True)
+
+# ----------------------------------------------------------------
+# setup color maps and opacity mapes used in the visualization
+# note: the Get..() functions create a new object, if needed
+# ----------------------------------------------------------------
+
+# get opacity transfer function/opacity map for 'RandomPointScalars'
+randomPointScalarsPWF = GetOpacityTransferFunction('RandomPointScalars')
+randomPointScalarsPWF.Points = [2.0, 0.0, 0.5, 0.0, 9.0, 1.0, 0.5, 0.0]
+randomPointScalarsPWF.ScalarRangeInitialized = 1
+
+# ----------------------------------------------------------------
+# finally, restore active source
+SetActiveSource(matlab_exportvtk)
+# ----------------------------------------------------------------
\ No newline at end of file

Matlab-Programs/TestCompute.m

+
+
+mu_1 = 10;
+rho_1 = 1;
+theta = 0.05;
+% alpha = -0.5;
+% alpha = 0.5;
+% alpha = 1;
+alpha = 18;
+beta = 0.5;  %egal welcher Wert Beta ... alha = 1 & Theta = 0.5 -> type 2 ? 
+beta = 2; 
+
+
+% Interessant: 
+alpha = 10;
+% alpha = 18;
+theta = 0.05;
+
+% unabhänhig von beta?
+beta = 2;
+
+% Compute components of B_eff
+% b1 = (mu_1*rho_1/4).*(beta./(theta+(1-theta).*beta)).*(1-theta.*(1+alpha));
+% b2 =  mu_1.*(rho_1/8).*(1-theta.*(1+beta.*alpha));
+
+
+%TEST (new)
+b1 = (3*rho_1/2).*beta.*(1-theta.*(1+alpha));
+b2 = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));
+
+
+
+mu_h = @(beta,theta) mu_1.*(beta./(theta+(1-theta).*beta));  % harmonic mean
+mu_bar = @(b,t) mu_1.*((1-theta)+theta.*beta);     % mu_bar
+
+
+
+%  q1 q2 q3..
+q1 = mu_h(beta,theta);
+q2 = mu_bar(beta,theta); 
+
+
+
+fprintf('--------------------')
+
+
+fprintf(' alpha:%d' , alpha)
+fprintf(' beta:%d ' , beta)
+fprintf(' theta:%d ', theta )
+fprintf('-------------------- \n')
+fprintf('q1*b1^2:')
+q1*b1^2 
+
+fprintf('q2*b2^2:')
+q2*b2^2
+
+
+fprintf('Test')
+(9/4)*(beta/(1+beta))*((1/2)-alpha+(1/2)*alpha^2) 
\ No newline at end of file

Matlab-Programs/TestWriteVTK.m

+
+
+x = 1:0.05:10;
+y = 1:0.05:10;
+z = 1:0.05:10;
+
+
+
+[X,Y,Z] = meshgrid(x,y,z);
+
+
+s = size(X,1);
+
+
+% A1 = repmat(1,size,3);
+% A2 = repmat(2,size,3);
+% A3 = rand(size,4)+ones(size,4);
+% % A3 = repmat(2,10,4);
+
+
+A1 = repmat(1,s,50);
+A2 = repmat(2,s,50);
+A3 = rand(s,81)+ones(s,81);
+% A3 = repmat(2,10,4);
+
+
+A  = [A1 A3 A2];
+B  = repmat(A,1,1,s);
+
+
+vtkwrite('wind.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type', B);
+
+
+
+% INPUT Parameters
+mu_1 = 1;
+rho_1 = 1;
+theta = 0.05; 
+alpha = 25; 
+beta= 20;
+
+set_mu_gamma = 'q1';
+% set_mu_gamma = 'q2';
+print_output= false;
+
+x = linspace(0,10,30);     %~alpha   % most interesting area for alpha >= 0 
+% y = linspace(0.05,100,121);     %~beta
+y = linspace(2,10,10);     %~beta %TEST
+y = 2;   % 2D Phase Diagram
+z = linspace(0.01,0.99,30);  %~theta              %theta = 0 , 1  -> parabolic case
+[X,Y,Z] = meshgrid(x,y,z);
+
+
+[A,angle_3D,V_3D,~] = arrayfun(@(a,b,theta)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X,Y,Z,'UniformOutput', false);
+
+
+angle_3D = cell2mat(angle_3D);
+V_3D = cell2mat(V_3D);
+
+vtkwrite('PhaseDiagram.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type',angle_3D);
+
+vtkwrite('PhaseDiagramType.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type',V_3D);
+
+
+
+% 2D Phase Diagram:
+% fix parameter (beta here) 
+% beta= 5;
+% 
+% x = linspace(-50,50,121);     %~alpha
+% z = linspace(0.01,0.99,121);  %~theta              %theta = 0 , 1  -> parabolic case
+% [X,Z] = meshgrid(x,z);
+% 
+% [A,angle_2D,V_2D,~] = arrayfun(@(a,theta)classifyMIN(mu_1,rho_1,a,beta,theta,set_mu_gamma,print_output),X,Z,'UniformOutput', false);
+% 
+% angle_2D = cell2mat(angle_2D);
+% V_2D = cell2mat(V_2D);
+% 
+
+% Y = ones(size(X))*beta;
+% vtkwrite('PhaseDiagram2D.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type',angle_2D);
+% 
+% vtkwrite('PhaseDiagramType2D.vtk', 'structured_grid', X,Y,Z, 'scalars', 'type',V_2D);
+
+
+
+

Matlab-Programs/classifyMIN.m

+
+function [A, angle, type, kappa] = classifyMIN (mu_1,rho_1,a,b,t,set_mu_gamma,print_output)
+
+% returns
+%   A : Matrix of basis coefficients [a1,a2,a3]
+%
+%   type : 
+%   Type of minimizer 1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
+%
+
+
+
+type = 0; % either 1,2,3,4 
+
+
+mu_h = @(b,t) mu_1.*(b./(t+(1-t).*b));  % harmonic mean
+mu_bar = @(b,t) mu_1.*((1-t)+t.*b);     % mu_bar
+
+if (set_mu_gamma == 'q1')
+  mu_gamma = @(b,t) mu_h(b,t);
+end
+if (set_mu_gamma == 'q2')
+  mu_gamma = @(b,t) mu_bar(b,t);
+end
+if (set_mu_gamma == 'm')
+  mu_gamma = @(b,t) 0.5*(mu_h(b,t) + mu_bar(b,t));
+end
+
+%  q1 q2 q3..
+q1 = mu_h(b,t);
+q2 = mu_bar(b,t); 
+q3 = mu_gamma(b,t);
+
+
+
+
+% values for q1,q2,q3 should be positiv
+% assert((q1 > 0 ) & (q2 > 0 ) & (q3 > 0), 'At least one of q1,q2 or q3 is not positive' )
+
+
+% Compute components of B_eff (old)
+% b1 = (mu_1*rho_1/4).*(b./(t+(1-t).*b)).*(1-t.*(1+a));
+% b2 =  mu_1.*(rho_1/8).*(1-t.*(1+b.*a));
+
+%TEST (new)
+b1 =  (3*rho_1/2).*b.*(1-t.*(1+a));
+b2 =  (3*rho_1./(4.*((1-t)+t.*b))).*(1-t.*(1+b.*a));
+
+
+% H = [q1 q3; q3 q2];         
+%check condition of H first
+% fprintf('condition number of Matrix H: %d \n', cond(H));
+
+
+CaseCount = 0; %check if more than one case ever happens
+
+epsilon = eps;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARABOLIC CASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% if abs(det(A)) < epsilon * min(abs(det(A)),0)  
+if abs(q1*q2-q3^2) < epsilon     
+    
+    fprintf('determinant equal zero (parabolic case)')
+    fprintf('SHOULD NOT HAPPEN')
+    
+   
+    
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ELLIPTIC CASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if (q1*q2-q3^2 > epsilon)
+%     fprintf('determinant greater than zero (elliptic case)');
+    
+%     a1_star =  (q2(b,t).*q1(b,t).*b1(a,b,t) - q3(b,t).*q2(b,t).*b2(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 
+%     a2_star =  (q1(b,t).*q2(b,t).*b2(a,b,t) - q3(b,t).*q1(b,t).*b1(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 
+    a1_star =  (q2.*q1.*b1 - q3.*q2.*b2)./(q1.*q2-q3.^2); 
+    a2_star =  (q1.*q2.*b2 - q3.*q1.*b1)./(q1.*q2-q3.^2);
+    
+    prod = a1_star*a2_star;
+    
+    
+    if(prod > epsilon) % (E1)  inside Lamba             %
+        % (a1_star,a2_star) is unique minimizer lies inside Lambda 
+        % therefore Minimizer not aligned with axes
+        
+%         fprintf('\n elliptic-case: (E1)');
+        a1 = a1_star;
+        a2 = a2_star;
+        type = 3;
+        CaseCount = CaseCount + 1;
+    end
+    % Make distinction between boundary & outside (prod < 0 ) 
+    if(abs(prod) < epsilon)  % (E2) on boundary of Lambda 
+%         fprintf('\n elliptic-case: (E2)');
+          
+        % Uniqueness of gloal minimizer if lies on boundary  if prod = 0 
+        % ----- % 
+        
+
+        % global minimizer lies on the boundary of Lambda depending on
+        % condition:
+        if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
+            a1 = b1;
+            a2 = 0*b1;
+            type = 1;  % Minimizer aligned with x1-axis
+            CaseCount = CaseCount + 1;
+        end
+        if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
+            a1 = 0*b1;
+            a2 = b2;
+            type = 2; % Minimizer aligned with x2-axis
+            CaseCount = CaseCount + 1;
+        end
+        if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
+            % two Minimizers ..pick one
+            a1 = b1;
+            a2 = 0*b1;
+            type = 4;
+            CaseCount = CaseCount + 1;
+        end
+    end
+    if((prod) < -1*epsilon) %Outside of Lambda 
+        
+        if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
+            a1 = b1;
+            a2 = 0*b1;
+            type = 1;  % Minimizer aligned with x1-axis
+            CaseCount = CaseCount + 1;
+        end
+        if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
+            a1 = 0*b1;
+            a2 = b2;
+            type = 2; % Minimizer aligned with x2-axis
+            CaseCount = CaseCount + 1;
+        end
+    end
+    
+end
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HYPERBOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if (q1*q2-q3^2 < -1*epsilon) 
+%     fprintf('determinant less than zero (hyperbolic case)');
+
+    if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
+        a1 = b1;
+        a2 = 0*b1;
+        type = 1;  % Minimizer aligned with x1-axis
+        CaseCount = CaseCount + 1;
+    end
+    if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
+        a1 = 0*b1;
+        a2 = b2;
+        type = 2; % Minimizer aligned with x2-axis
+        CaseCount = CaseCount + 1;
+    end
+    if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
+        % two Minimizers ..pick one
+        a1 = b1;
+        a2 = 0*b1;
+        type = 4;
+        CaseCount = CaseCount + 1;
+    end
+
+    
+    % CAN NOT BE TYPE 3!!
+    
+end
+
+
+
+% Compute a3 from a1 % a2
+a3 = sqrt(2*a1*a2);                    % WRONG ? 
+
+
+
+
+% compute angle between [sqrt(a1) , sqrt(a2)] and e1:
+% angle = atan2(sqrt(a2),sqrt(a1));
+
+v = [sqrt(a1), sqrt(a2)];
+e1 = [1 0];
+
+if (type == 3 )
+%    angle = atan2(a2,a1);
+%    angle = acos(v*e1/(sqrt(a1+a2)));
+else
+   angle = 0;
+end
+% angle = atan2(norm(cross(a,b)), dot(a,b))
+
+% Try to compute angle of Matrices
+
+E1 = [1 0;0 0 ];
+% V = [a1 sqrt(a1*a2)/sqrt(2); sqrt(a1*a2)/sqrt(2) a2]; % This ?? NO
+
+V = [a1 sqrt(a1*a2); sqrt(a1*a2) a2];         % This!
+
+% CHECK
+% U = trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1)))   
+% if abs(U) > 1
+%    fprintf('value greater 1') 
+% end
+
+
+% angle = atan2(sqrt(abs(a2)), sqrt(abs(a1)));   % does this make sense? 
+
+
+
+% angle = acos(trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1))));    %angle in radians
+% angle = acos(trace(V'*E1)/(sqrt(trace(V'*V))*sqrt(trace(E1'*E1))))/pi * 180; % angle in degrees
+
+
+% CHECK does case (0,-b) ever happen?  Yes 
+% if  (q2*b2^2 > q1*b1^2 && b2 < 0)
+%     fprintf('point lies on negative half of y-axis'); 
+% end
+
+
+
+% CHECK if Minimizer ever lies inside lambda ... YES 
+% if  (a1 ~= 0 && a2 ~= 0)
+%     fprintf('minimizer lies inside lambda'); 
+% end
+
+
+% Alternative atan2d(y,x): returns angle in degrees 
+
+
+% CHECK
+e = [sqrt(a1), sqrt(a2)];        % Might be complex...
+norm = sqrt(a1+a2);              % Might be complex...
+e = e./norm;
+angle = atan2(e(2), e(1));   %TEST
+
+if (imag(e(1))~= 0 || imag(e(2))~=0)
+   fprintf('complex vector e'); 
+end
+
+% if (imag(norm)~=0)
+%    fprintf('complex norm'); % happens a lot .. 
+% end
+
+
+% do it this way to avoid complex numbers: 
+e = [sqrt((a1/(a1+a2))), sqrt((a2/(a1+a2)))];    % always positive under sqrt here .. basically takes absolute value here
+angle = atan2(e(2), e(1)); 
+
+
+if (imag(e(1))~= 0 || imag(e(2))~=0)
+   fprintf('complex vector e'); 
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TEST  -------------------
+% the magnitude kappa of the vector does not matter for the angle ..
+% In order to also get negative angles just use:
+% atan2(a2,a1) 
+
+% angle = atan2(a2,a1);         %% FEHLER : muss [sqrt(a1) sqrt(a2)] betrachten..
+% angle = atan2(sqrt(a2),sqrt(a1));  % Inputs must be real... 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+% angle2 = acos(e(1)/(sqrt(e(1)^2 + e(2)^2) ));
+% if angle ~= angle2
+%    interesting = 1; 
+% end
+
+% if (angle ~= 0 | angle ~= 3.1416)
+%    fprintf('angle not zero or pi.. plot type')
+%    type  
+%    angle
+% end
+
+%compute Kappa? 
+% fprintf('Output Kappa:')
+% k = sqrt(abs(a1) + abs(a2));  % ? 
+kappa = (a1 + a2); 
+
+
+
+% e = [sqrt(a1), sqrt(a2)];
+% norm = sqrt((a1+a2));
+% e = e./norm;
+% K.*(e'*e);
+
+
+if (CaseCount > 1)
+  fprintf('=============================== ERROR========================================= \n')
+  fprintf('Error: More than one Case happened!')
+end
+
+
+
+% Coefficients of minimizer 
+
+if(print_output)
+   fprintf(' \n ') 
+   fprintf('=============================== OUTPUT ========================================= \n')
+   fprintf('Minimizer is of Type: %d  \n' , type);
+   fprintf('Coefficients a1,a2,a3 given by : %d, %d, %d   \n', a1, a2, a3);
+   fprintf('================================================================================ \n')
+end
+
+
+A = [a1, a2, a3];
+
+end

Matlab-Programs/compute_F.m

+function F = compute_F(alpha,B,q1,q2,q3)
+
+% r = compute_r(alpha,B,q1,q2,q3);
+
+v = [cos(alpha);sin(alpha)];
+
+b1 = B(1,1);
+b2 = B(2,2);
+b3 = B(1,2);
+
+%compute Q(v_alpha x v_alpha)
+Q = q1.*v(1).^4 + q2.*v(2).^4 + 2.*q3.*v(1).^2.*v(2).^2;
+% 
+% TP = v*v';
+% L = stVenant(TP,mu,lambda);
+
+tmp1 = q1.*(v(1).^2+b1).^2 + q2.*(v(2).^2+b2).^2 + q3.*(sqrt(2)*v(1)*v(2)+b3).^2;
+tmp2 = q1.*b1.^2 + q2.*b2.^2+ q3.*b3.^2;
+L = 0.5*(tmp1-Q-tmp2) ;       %Polarization identity 
+
+
+r = L./Q;
+
+
+% F = r.^2.*Q - 2.*r.*trace(L'*B)
+F = (r.^2).*Q - 2.*r.*L;
+
+
+end
+

Matlab-Programs/compute_r.m

+function r = compute_r(alpha,B,q1,q2,q3)
+
+% r = compute_r(alpha,B,q1,q2,q3);
+
+v = [cos(alpha);sin(alpha)];
+
+b1 = B(1,1);
+b2 = B(2,2);
+b3 = B(1,2);
+
+%compute Q(v_alpha x v_alpha)
+Q = q1.*v(1).^4 + q2.*v(2).^4 + 2.*q3.*v(1).^2.*v(2).^2;
+
+
+tmp1 = q1.*(v(1).^2+b1).^2 + q2.*(v(2).^2+b2).^2 + q3.*(sqrt(2)*v(1)*v(2)+b3).^2;
+tmp2 = q1.*b1.^2 + q2.*b2.^2+ q3.*b3.^2;
+L = 0.5*(tmp1-Q-tmp2) ;       %Polarization identity 
+
+r = L./Q;
+
+
+
+
+end
+

Matlab-Programs/findParabolic.m

+
+
+syms f(alpha,beta,theta)
+syms b1(alpha,beta,theta)
+syms b2(alpha,beta,theta)
+syms q1(beta,theta)
+syms q2(beta,theta)
+
+mu_1 = sym(1);
+rho_1 = sym(1);
+
+%TEST (new)
+b1(alpha,beta,theta) = (3*rho_1/2).*beta.*(1-theta.*(1+alpha));
+b2(alpha,beta,theta) = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));
+
+q1(beta,theta) = mu_1.*(beta./(theta+(1-theta).*beta)) ;
+q2(beta,theta) = mu_1.*((1-theta)+theta.*beta); 
+
+
+
+fprintf('Functions to be minimized')
+% f(beta,theta,alpha) = (mu_1.*(beta./(theta+(1-theta).*beta))).*((3*rho_1/2).*beta.*(1-theta.*(1+alpha))).^2 - (mu_1.*((1-theta)+theta.*beta)).*(3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha)).^2;
+f(beta,theta,alpha) = q1(beta,theta).*b1(alpha,beta,theta).^2 - (q2(beta,theta).*b2(alpha,beta,theta).^2) ;
+
+% Fix Theta ... 
+thetav = sym(1/2);
+f = subs(f,theta,thetav);
+b1 = subs(b1,theta,thetav);
+b2 = subs(b2,theta,thetav);
+q1 = subs(q1,theta,thetav);
+q2 = subs(q2,theta,thetav);
+
+
+
+% f = subs(f,theta,0.25);
+% b1 = subs(b1,theta,0.25);
+% b2 = subs(b2,theta,0.25);
+% q1 = subs(q1,theta,0.25);
+% q2 = subs(q2,theta,0.25);
+
+
+
+% f = subs(f,alpha,3);
+
+
+eq = f == 0;
+
+
+S = solve(eq,alpha, beta,'MaxDegree' , 5)
+% S = solve(eq,alpha, beta,theta,'MaxDegree' , 5)
+
+S.alpha
+S.beta
+
+
+% CHECK
+b1 = subs(b1,alpha,S.alpha(1));
+b1 = subs(b1,beta,S.beta(1));
+b2 = subs(b2,alpha,S.alpha(1));
+b2 = subs(b2,beta,S.beta(1));
+
+q1 = subs(q1,beta,S.beta(1));
+q2 = subs(q2,beta,S.beta(1));
+
+b1 = double(b1)
+b2 = double(b2)
+q1 = double(q1)
+q2 = double(q2)
+
+
+b1.*((q1).^2)- ( b2.*(q2.^2)) 
+
+
+
+f = subs(f,alpha,S.alpha(1));
+f = subs(f,beta,S.beta(1))
+simplify(f)
+
+
+% 
+% 
+% mu_h = @(beta,theta) mu_1.*(beta./(theta+(1-theta).*beta));  % harmonic mean
+% mu_bar = @(b,t) mu_1.*((1-theta)+theta.*beta);     % mu_bar
+% 
+% 
+% 
+% %  q1 q2 q3..
+% q1 = mu_h(beta,theta);
+% q2 = mu_bar(beta,theta); 
+
+
+
+fprintf('--------------------')
+
+

Matlab-Programs/findParameters.m

+
+
+
+syms alpha beta theta rho_1 mu_1 
+assume(beta,'positive')
+assume(theta,'positive')
+assume(mu_1,'positive')
+assume(rho_1,'positive')
+
+% syms f(alpha,beta,theta)
+% syms b1(alpha,beta,theta)
+% syms b2(alpha,beta,theta)
+% syms q1(beta,theta)
+% syms q2(beta,theta)
+
+% mu_1 = sym(1);
+% rho_1 = sym(1);
+
+%TEST (new)
+b1(alpha,beta,theta,rho_1) = (3*rho_1/2).*beta.*(1-theta.*(1+alpha));
+b2(alpha,beta,theta,rho_1) = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));
+
+
+
+
+q1(beta,theta,mu_1) = mu_1.*(beta./(theta+(1-theta).*beta)) ;
+q2(beta,theta,mu_1) = mu_1.*((1-theta)+theta.*beta); 
+q3_star(beta,theta,mu_1) = (q1(beta,theta,mu_1)*q2(beta,theta,mu_1))^(1/2)
+
+% 
+% fprintf('Functions to be minimized')
+% % f(beta,theta,alpha) = (mu_1.*(beta./(theta+(1-theta).*beta))).*((3*rho_1/2).*beta.*(1-theta.*(1+alpha))).^2 - (mu_1.*((1-theta)+theta.*beta)).*(3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha)).^2;
+% f(beta,theta,alpha) = q1(beta,theta).*b1(alpha,beta,theta).^2 - (q2(beta,theta).*b2(alpha,beta,theta).^2) ;
+
+% Fix Theta ... 
+% thetav = sym(1/4);
+% q3_star = subs(q3_star,theta,thetav);
+% b1 = subs(b1,theta,thetav);
+% b2 = subs(b2,theta,thetav);
+% q1 = subs(q1,theta,thetav);
+% q2 = subs(q2,theta,thetav);
+
+
+% theta = thetav
+eq1 = (3*rho_1/2).*beta.*(1-theta.*(1+alpha)) == (q1(beta,theta,mu_1)*q2(beta,theta,mu_1))^(1/2)
+
+eq2 = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha)) == mu_1.*(beta./(theta+(1-theta).*beta)) 
+
+eqns = [eq1, eq2] 
+
+% eqns = subs(eqns,{theta,alpha,rho_1},{1/20,2,1} ) % works! 
+% eqns = subs(eqns,{rho_1,alpha,theta},{1,2,1/4} )
+eqns = subs(eqns,{theta},{1/4} )                        % fix only theta...
+% eqns = subs(eqns,{theta,mu_1,rho_1},{1/4,10,1} )
+
+
+% eqns = subs(eqns,{mu_1,beta},{387/1250,7741/10000} )
+% simplify(eqns)
+
+% eqns = subs(eqns,{theta,mu_1,rho_1},{1/4,10,5} )
+% eqns = subs(eqns,beta,2.0)
+
+% S = solve(eqns,mu_1,beta,'MaxDegree' , 5)
+% S = solve(eqns,theta,beta,'MaxDegree' , 5)
+% S = solve(eqns,theta,mu_1,'MaxDegree' , 5)
+
+
+%interesting:
+% eqns = subs(eqns,{theta,alpha,rho_1},{1/20,2,1} )
+% S = solve(eqns,mu_1,beta,'MaxDegree' , 5)
+
+
+% S.mu_1
+% S.beta
+% double(S.mu_1)
+% double(S.beta)
+% 
+% 
+% eqns = subs(eqns,{mu_1,beta},{0.3096,0.7741} )
+% simplify(eqns)
+
+
+S = solve(eqns, alpha,beta, rho_1,mu_1, 'MaxDegree' , 5)
+% S = solve(eqns, mu_1,beta, 'MaxDegree' , 5)
+
+% S = solve(eqns, alpha, beta, 'MaxDegree' , 5)
+
+S.alpha
+S.beta
+double(S.alpha)
+double(S.beta)
+
+eqns = subs(eqns,{alpha,beta},{95.2970,0.1577} )
+% simplify(eqns)
+
+
+% f = subs(f,theta,0.25);
+% b1 = subs(b1,theta,0.25);
+% b2 = subs(b2,theta,0.25);
+% q1 = subs(q1,theta,0.25);
+% q2 = subs(q2,theta,0.25);
+
+
+
+% f = subs(f,alpha,3);
+
+
+
+
+% S = solve(eq,alpha, beta,'MaxDegree' , 5)
+% S = solve(eq,alpha, beta,theta,'MaxDegree' , 5)
+
+% % S.alpha
+% S.beta
+% % 
+% 
+% % CHECK
+% b1 = subs(b1,alpha,S.alpha(1));
+% b1 = subs(b1,beta,S.beta(1));
+% b2 = subs(b2,alpha,S.alpha(1));
+% b2 = subs(b2,beta,S.beta(1));
+% 
+% q1 = subs(q1,beta,S.beta(1));
+% q2 = subs(q2,beta,S.beta(1));
+% 
+% b1 = double(b1)
+% b2 = double(b2)
+% q1 = double(q1)
+% q2 = double(q2)
+% 
+% 
+% b1.*((q1).^2)- ( b2.*(q2.^2)) 
+% 
+% 
+% 
+% f = subs(f,alpha,S.alpha(1));
+% f = subs(f,beta,S.beta(1))
+% simplify(f)
+
+
+% 
+% 
+% mu_h = @(beta,theta) mu_1.*(beta./(theta+(1-theta).*beta));  % harmonic mean
+% mu_bar = @(b,t) mu_1.*((1-theta)+theta.*beta);     % mu_bar
+% 
+% 
+% 
+% %  q1 q2 q3..
+% q1 = mu_h(beta,theta);
+% q2 = mu_bar(beta,theta); 
+
+
+
+fprintf('--------------------')
+
+

Matlab-Programs/matlab.txt

+1 1 1.906127239467446
+1 2 -2.84445771411170261e-27
+1 3 1.36674531873130088e-11
+2 1 -1.70707845907209275e-26
+2 2 2.0833333333333246
+2 3 6.44949347783633458e-19
+3 1 2.21186925202339457e-10
+3 2 6.4495181619832732e-19
+3 3 1.92304720901783655

Matlab-Programs/matlab_export.vtk

+# vtk DataFile Version 2.0
+VTK from Matlab
+ASCII
+DATASET POLYDATA
+POINTS 102 float
+1.00 0.84 1.00 2.00 0.91 1.41 3.00 0.14 1.73 
+4.00 -0.76 2.00 5.00 -0.96 2.24 6.00 -0.28 2.45 
+7.00 0.66 2.65 8.00 0.99 2.83 9.00 0.41 3.00 
+10.00 -0.54 3.16 11.00 -1.00 3.32 12.00 -0.54 3.46 
+13.00 0.42 3.61 14.00 0.99 3.74 15.00 0.65 3.87 
+16.00 -0.29 4.00 17.00 -0.96 4.12 18.00 -0.75 4.24 
+19.00 0.15 4.36 20.00 0.91 4.47 21.00 0.84 4.58 
+22.00 -0.01 4.69 23.00 -0.85 4.80 24.00 -0.91 4.90 
+25.00 -0.13 5.00 26.00 0.76 5.10 27.00 0.96 5.20 
+28.00 0.27 5.29 29.00 -0.66 5.39 30.00 -0.99 5.48 
+31.00 -0.40 5.57 32.00 0.55 5.66 33.00 1.00 5.74 
+34.00 0.53 5.83 35.00 -0.43 5.92 36.00 -0.99 6.00 
+37.00 -0.64 6.08 38.00 0.30 6.16 39.00 0.96 6.24 
+40.00 0.75 6.32 41.00 -0.16 6.40 42.00 -0.92 6.48 
+43.00 -0.83 6.56 44.00 0.02 6.63 45.00 0.85 6.71 
+46.00 0.90 6.78 47.00 0.12 6.86 48.00 -0.77 6.93 
+49.00 -0.95 7.00 50.00 -0.26 7.07 51.00 0.67 7.14 
+52.00 0.99 7.21 53.00 0.40 7.28 54.00 -0.56 7.35 
+55.00 -1.00 7.42 56.00 -0.52 7.48 57.00 0.44 7.55 
+58.00 0.99 7.62 59.00 0.64 7.68 60.00 -0.30 7.75 
+61.00 -0.97 7.81 62.00 -0.74 7.87 63.00 0.17 7.94 
+64.00 0.92 8.00 65.00 0.83 8.06 66.00 -0.03 8.12 
+67.00 -0.86 8.19 68.00 -0.90 8.25 69.00 -0.11 8.31 
+70.00 0.77 8.37 71.00 0.95 8.43 72.00 0.25 8.49 
+73.00 -0.68 8.54 74.00 -0.99 8.60 75.00 -0.39 8.66 
+76.00 0.57 8.72 77.00 1.00 8.77 78.00 0.51 8.83 
+79.00 -0.44 8.89 80.00 -0.99 8.94 81.00 -0.63 9.00 
+82.00 0.31 9.06 83.00 0.97 9.11 84.00 0.73 9.17 
+85.00 -0.18 9.22 86.00 -0.92 9.27 87.00 -0.82 9.33 
+88.00 0.04 9.38 89.00 0.86 9.43 90.00 0.89 9.49 
+91.00 0.11 9.54 92.00 -0.78 9.59 93.00 -0.95 9.64 
+94.00 -0.25 9.70 95.00 0.68 9.75 96.00 0.98 9.80 
+97.00 0.38 9.85 98.00 -0.57 9.90 99.00 -1.00 9.95 
+100.00 -0.51 10.00 0.00 0.00 0.00 0.00 0.00 0.00 
+
+LINES 198 594
+2 0 1
+2 1 2
+2 2 3
+2 3 4
+2 4 5
+2 5 6
+2 6 7
+2 7 8
+2 8 9
+2 9 10
+2 10 11
+2 11 12
+2 12 13
+2 13 14
+2 14 15
+2 15 16
+2 16 17
+2 17 18
+2 18 19
+2 19 20
+2 20 21
+2 21 22
+2 22 23
+2 23 24
+2 24 25
+2 25 26
+2 26 27
+2 27 28
+2 28 29
+2 29 30
+2 30 31
+2 31 32
+2 32 33
+2 33 34
+2 34 35
+2 35 36
+2 36 37
+2 37 38
+2 38 39
+2 39 40
+2 40 41
+2 41 42
+2 42 43
+2 43 44
+2 44 45
+2 45 46
+2 46 47
+2 47 48
+2 48 49
+2 49 50
+2 50 51
+2 51 52
+2 52 53
+2 53 54
+2 54 55
+2 55 56
+2 56 57
+2 57 58
+2 58 59
+2 59 60
+2 60 61
+2 61 62
+2 62 63
+2 63 64
+2 64 65
+2 65 66
+2 66 67
+2 67 68
+2 68 69
+2 69 70
+2 70 71
+2 71 72
+2 72 73
+2 73 74
+2 74 75
+2 75 76
+2 76 77
+2 77 78
+2 78 79
+2 79 80
+2 80 81
+2 81 82
+2 82 83
+2 83 84
+2 84 85
+2 85 86
+2 86 87
+2 87 88
+2 88 89
+2 89 90
+2 90 91
+2 91 92
+2 92 93
+2 93 94
+2 94 95
+2 95 96
+2 96 97
+2 97 98
+2 98 99
+2 1 0
+2 2 1
+2 3 2
+2 4 3
+2 5 4
+2 6 5
+2 7 6
+2 8 7
+2 9 8
+2 10 9
+2 11 10
+2 12 11
+2 13 12
+2 14 13
+2 15 14
+2 16 15
+2 17 16
+2 18 17
+2 19 18
+2 20 19
+2 21 20
+2 22 21
+2 23 22
+2 24 23
+2 25 24
+2 26 25
+2 27 26
+2 28 27
+2 29 28
+2 30 29
+2 31 30
+2 32 31
+2 33 32
+2 34 33
+2 35 34
+2 36 35
+2 37 36
+2 38 37
+2 39 38
+2 40 39
+2 41 40
+2 42 41
+2 43 42
+2 44 43
+2 45 44
+2 46 45
+2 47 46
+2 48 47
+2 49 48
+2 50 49
+2 51 50
+2 52 51
+2 53 52
+2 54 53
+2 55 54
+2 56 55
+2 57 56
+2 58 57
+2 59 58
+2 60 59
+2 61 60
+2 62 61
+2 63 62
+2 64 63
+2 65 64
+2 66 65
+2 67 66
+2 68 67
+2 69 68
+2 70 69
+2 71 70
+2 72 71
+2 73 72
+2 74 73
+2 75 74
+2 76 75
+2 77 76
+2 78 77
+2 79 78
+2 80 79
+2 81 80
+2 82 81
+2 83 82
+2 84 83
+2 85 84
+2 86 85
+2 87 86
+2 88 87
+2 89 88
+2 90 89
+2 91 90
+2 92 91
+2 93 92
+2 94 93
+2 95 94
+2 96 95
+2 97 96
+2 98 97
+2 99 98

Matlab-Programs/redblue.m

+function c = redblue(m)
+%REDBLUE    Shades of red and blue color map
+%   REDBLUE(M), is an M-by-3 matrix that defines a colormap.
+%   The colors begin with bright blue, range through shades of
+%   blue to white, and then through shades of red to bright red.
+%   REDBLUE, by itself, is the same length as the current figure's
+%   colormap. If no figure exists, MATLAB creates one.
+%
+%   For example, to reset the colormap of the current figure:
+%
+%             colormap(redblue)
+%
+%   See also HSV, GRAY, HOT, BONE, COPPER, PINK, FLAG, 
+%   COLORMAP, RGBPLOT.
+
+%   Adam Auton, 9th October 2009
+
+if nargin < 1, m = size(get(gcf,'colormap'),1); end
+
+if (mod(m,2) == 0)
+    % From [0 0 1] to [1 1 1], then [1 1 1] to [1 0 0];
+    m1 = m*0.5;
+    r = (0:m1-1)'/max(m1-1,1);
+    g = r;
+    r = [r; ones(m1,1)];
+    g = [g; flipud(g)];
+    b = flipud(r);
+else
+    % From [0 0 1] to [1 1 1] to [1 0 0];
+    m1 = floor(m*0.5);
+    r = (0:m1-1)'/max(m1,1);
+    g = r;
+    r = [r; ones(m1+1,1)];
+    g = [g; 1; flipud(g)];
+    b = flipud(r);
+end
+
+c = [r g b]; 
+

Matlab-Programs/resources/addons_core.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<addonCore>
+  <label>Red Blue Colormap</label>
+  <version>1.0.0.0</version>
+  <type>zip</type>
+  <identifier>e5698820-4a80-11e4-9553-005056977bd0</identifier>
+  <createdBy name="Adam Auton"/>
+  <image>resources/screenshot.png</image>
+</addonCore>

Matlab-Programs/resources/matlab_path_entries.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<paths>
+  <path>.</path>
+</paths>

Matlab-Programs/resources/metadata.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<addOn>
+  <identifier>e5698820-4a80-11e4-9553-005056977bd0</identifier>
+  <displayType>Function</displayType>
+  <translatedDisplayType>
+    <en_US>Function</en_US>
+    <ja_JP>関数</ja_JP>
+    <ko_KR>함수</ko_KR>
+    <zh_CN>函数</zh_CN>
+  </translatedDisplayType>
+  <name>Red Blue Colormap</name>
+  <author>Adam Auton</author>
+  <version>1.0.0.0</version>
+  <downloadUrl>https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/25536/versions/1/download/zip?src=addons_ml_desktop_install&amp;profile_id=14169257&amp;license=40758619&amp;release_family=R2020b</downloadUrl>
+  <licenseUrl>https://addons.mathworks.com/registry/v1/e5698820-4a80-11e4-9553-005056977bd0/1.0.0.0/-/license</licenseUrl>
+  <previewImageUrl>https://www.mathworks.com/responsive_image/160/120/0/0/0/cache/matlabcentral/mlc-downloads/downloads/submissions/25536/versions/1/screenshot.png</previewImageUrl>
+  <releaseNotesUrl>https://www.mathworks.com/add-ons/e5698820-4a80-11e4-9553-005056977bd0/d1034754-ea44-08f8-b658-22b590dfae7e/releaseNotes</releaseNotesUrl>
+  <installationFolder>Functions</installationFolder>
+</addOn>