PolyGridWizZ/Poly6_end.m at master · ZaunerM/PolyGridWizZ · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
function s_C1 = Poly6_end(N,s1,s2,ds1,ds2,dds1,dds2,ddds2,plotting)
bx=[1:N];
%bx=linspace(x1,x2,N);
dx=1; %deriv(bx',1);

syms cc dd ee ff gg hh ii
%%x
eqn1 =       cc * bx(1)^6 +      dd * bx(1)^5 +      ee * bx(1)^4 +      ff * bx(1)^3 +    gg * bx(1)^2 +   hh * bx(1)^1 +ii == s1;
eqn2 =       cc*bx(end)^6 +      dd*bx(end)^5 +      ee*bx(end)^4 +      ff*bx(end)^3 +    gg*bx(end)^2 +   hh*bx(end)^1 +ii == s2;
%%spacing
eqn3 =     6*cc * bx(1)^5 +    5*dd * bx(1)^4 +    4*ee * bx(1)^3 +    3*ff * bx(1)^2 +  2*gg * bx(1)^1 +1*hh           +0  == ds1/dx(1);
eqn4 =     6*cc*bx(end)^5 +    5*dd*bx(end)^4 +    4*ee*bx(end)^3 +    3*ff*bx(end)^2 +  2*gg*bx(end)^1 +1*hh           +0  == ds2/dx(end);
%%second derivative of x
eqn7 =   5*6*cc * bx(1)^4 +  4*5*dd * bx(1)^3 +  3*4*ee * bx(1)^2 +  2*3*ff * bx(1)^1 +1*2*gg           +0              +0  == dds1/(dx(1)^2);
eqn5 =   5*6*cc*bx(end)^4 +  4*5*dd*bx(end)^3 +  3*4*ee*bx(end)^2 +  2*3*ff*bx(end)^1 +1*2*gg           +0              +0  == dds2/(dx(end)^2);
%%3rd deriv
eqn6 = 4*5*6*cc * bx(end)^3 +3*4*5*dd * bx(end)^2 +2*3*4*ee * bx(end)^1 +1*2*3*ff           +0                +0              +0  == ddds2/dx(end);

[As,Bs] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7], [cc,dd,ee,ff,gg,hh,ii]);

%Sol1 = linsolve(As,Bs);
Sol1 = mldivide(As,Bs);

s_C1_=double(Sol1(1,1))*bx.^6+double(Sol1(2,1))*bx.^5+double(Sol1(3,1))*bx.^4+double(Sol1(4,1))*bx.^3+double(Sol1(5,1))*bx.^2+double(Sol1(6,1))*bx.^1+double(Sol1(7,1))*bx.^0;
s_C1=s_C1_;

if plotting=='t'
    figure
    %plot(deriv(deriv(s_C1',1)./dx,1)./dx)
    hold on
    plot(deriv(s_C1_',1),'r')
%     plot(s_C1_','r')
end

% r1=linspace(0,(N-1)*ds1,N);
% r2=linspace(s2,s2-(N-1)*ds2,N);
%
% lam=[1:Nb];
% a=-2/( 2*((1-Nb^3)-3*(1-Nb)) - ((1-Nb^2)-2*(1-Nb))*3*(1+Nb) );
% b=-a*(3*(1+Nb))/2;
% c=-3*a-2*b;
% d=-a-b-c;
% bl_=a*lam.^3+b*lam.^2+c*lam+d;
% bl(1)=0;
% bl(2:Nb+1)=bl_;
% bl(Nb+2)=1;
% s_C1=s_C1_;
% for i=1:Nb+2
%     s_C1(i)=r1(i)+bl(i)*(s_C1_(i)-r1(i));
%     s_C1(end+1-i)=r2(i)+bl(i)*(s_C1_(end+1-i)-r2(i));
% end
% if plotting=='t'
%    plot(deriv(s_C1',1),'k')
% %        plot(s_C1','k')
% %        plot(r1','b')
% end
% figure
% plot(r1)
% hold on
% plot(s_C1_)
% plot(s_C1)

end