Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger1dCurve.m @ 511:57f3493f851b feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added sqrt of Ji in the right places, not sure about the interfaces, will not test it properly now
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 08 Jun 2017 10:33:36 +0200
parents 508b7493be94
children 32a24485f3e8
comparison
equal deleted inserted replaced
510:71908bbce2e1 511:57f3493f851b
84 % neighbour_scheme is an instance of Scheme that should be interfaced to. 84 % neighbour_scheme is an instance of Scheme that should be interfaced to.
85 % neighbour_boundary is a string specifying which boundary to interface to. 85 % neighbour_boundary is a string specifying which boundary to interface to.
86 86
87 function [D] = d_fun(obj,t) 87 function [D] = d_fun(obj,t)
88 obj.variable_update(t); % In driscretization? 88 obj.variable_update(t); % In driscretization?
89 D = (-0.5*(obj.D1*obj.a + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat); 89 D = sqrt(obj.Ji)*(-0.5*(obj.D1*obj.a + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat)*sqrt(obj.Ji);
90 % D = (-0.5*(obj.D1*obj.a -obj.a_xi+ obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat); 90 % D = (-0.5*(obj.D1*obj.a -obj.a_xi+ obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
91 % D= obj.Ji*(-sqrt(obj.a)*obj.D1*sqrt(obj.a) + 0.5*obj.a_xi + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat); 91 % D= obj.Ji*(-sqrt(obj.a)*obj.D1*sqrt(obj.a) + 0.5*obj.a_xi + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
92 end 92 end
93 93
94 94
119 switch type 119 switch type
120 % Dirichlet boundary condition 120 % Dirichlet boundary condition
121 case {'D','d','dirichlet'} 121 case {'D','d','dirichlet'}
122 tau1 = @(t) s * 1i*obj.Ji(p,p)*d; 122 tau1 = @(t) s * 1i*obj.Ji(p,p)*d;
123 tau2 = @(t) (1*s*obj.a(p,p))/2*e; 123 tau2 = @(t) (1*s*obj.a(p,p))/2*e;
124 closure = @(t)(obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e'); 124 closure = @(t)obj.Hi*sqrt(obj.Ji)*(tau1(t)*e' + tau2(obj.a)*e')*sqrt(obj.Ji);
125 125
126 switch class(data) 126 switch class(data)
127 case 'double' 127 case 'double'
128 penalty = @(t) -(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data); 128 penalty = @(t) -obj.Hi*sqrt(obj.Ji)*(tau1*data+tau2(obj.a)*data)*sqrt(obj.Ji);
129 % case 'function_handle' 129 % case 'function_handle'
130 % penalty = @(t)-obj.Hi*tau*data(t); 130 % penalty = @(t)-obj.Hi*tau*data(t);
131 otherwise 131 otherwise
132 error('Weird data argument!') 132 error('Weird data argument!')
133 end 133 end
144 [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary); 144 [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary);
145 [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 145 [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
146 146
147 a1 = s_u* 1/2 * 1i ; 147 a1 = s_u* 1/2 * 1i ;
148 b1 = -s_u* 1/2 * 1i; 148 b1 = -s_u* 1/2 * 1i;
149 gamma = @(a) -s_u*a(p_u,p_u)/2*e_u; 149 gamma = @(a) s_u*a(p_u,p_u)/2*e_u;
150 150
151 tau = @(t) a1*obj.Ji(p_u,p_u)^2*d_u; 151 tau = @(t) a1*obj.Ji(p_u,p_u)^2*d_u;
152 sig = b1*e_u; 152 sig = b1*e_u;
153 153
154 closure = @(t) obj.Hi * (tau(t)*e_u' + sig*obj.Ji(p_u,p_u)^2*d_u' + gamma(obj.a)*e_u'); 154 closure = @(t) obj.Hi * (tau(t)*e_u' + sig*obj.Ji(p_u,p_u)^2*d_u' + obj.Ji(p_u,p_u)*gamma(obj.a)*e_u');
155 penalty = @(t) obj.Hi * (-tau(t)*e_v' - sig*obj.Ji(p_u,p_u)^2*d_v' - gamma(obj.a)*e_v'); 155 penalty = @(t) obj.Hi * (-tau(t)*e_v' - sig*obj.Ji(p_u,p_u)^2*d_v' - obj.Ji(p_u,p_u)*gamma(obj.a)*e_v');
156 end 156 end
157 157
158 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 158 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
159 % The right boundary is considered the positive boundary 159 % The right boundary is considered the positive boundary
160 function [e,d,s,p] = get_boundary_ops(obj,boundary) 160 function [e,d,s,p] = get_boundary_ops(obj,boundary)