Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_4.m @ 427:a613960a157b feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged with feature/beams
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 26 Jan 2017 15:59:25 +0100
parents 8b8672134be8
children
comparison
equal deleted inserted replaced
426:29944ea7674b 427:a613960a157b
1 function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_4(m,h) 1 function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_4(m,h)
2 2
3 BP = 4; 3 BP = 4;
4 if(m<2*BP) 4 if(m<2*BP)
5 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 5 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
6 end 6 end
7 7
18 HI=inv(H); 18 HI=inv(H);
19 19
20 20
21 Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12; 21 Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12;
22 %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); 22 %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));
23 Q_U = [0 0.59e2 / 0.96e2 -0.1e1 / 0.12e2 -0.1e1 / 0.32e2; -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2; 0.1e1 / 0.32e2 0 -0.59e2 / 0.96e2 0;]; 23 Q_U = [
24 0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2;
25 -0.59e2/0.96e2 0 0.59e2/0.96e2 0;
26 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2;
27 0.1e1/0.32e2 0 -0.59e2/0.96e2 0;
28 ];
24 Q(1:4,1:4)=Q_U; 29 Q(1:4,1:4)=Q_U;
25 Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 ); 30 Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 );
26 31
27 D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; 32 D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
28 33
29 M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12; 34 M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12;
30 35
31 %M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); 36 %M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0));
32 37
33 M_U=[0.9e1 / 0.8e1 -0.59e2 / 0.48e2 0.1e1 / 0.12e2 0.1e1 / 0.48e2; -0.59e2 / 0.48e2 0.59e2 / 0.24e2 -0.59e2 / 0.48e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.48e2 0.55e2 / 0.24e2 -0.59e2 / 0.48e2; 0.1e1 / 0.48e2 0 -0.59e2 / 0.48e2 0.59e2 / 0.24e2;]; 38 M_U=[
39 0.9e1/0.8e1 -0.59e2/0.48e2 0.1e1/0.12e2 0.1e1/0.48e2;
40 -0.59e2/0.48e2 0.59e2/0.24e2 -0.59e2/0.48e2 0;
41 0.1e1/0.12e2 -0.59e2/0.48e2 0.55e2/0.24e2 -0.59e2/0.48e2;
42 0.1e1/0.48e2 0 -0.59e2/0.48e2 0.59e2/0.24e2;
43 ];
34 M(1:4,1:4)=M_U; 44 M(1:4,1:4)=M_U;
35 45
36 M(m-3:m,m-3:m)=rot90( M_U ,2 ); 46 M(m-3:m,m-3:m)=rot90( M_U ,2 );
37 M=M/h; 47 M=M/h;
38 48
39 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h; 49 S_U=[-0.11e2/0.6e1 3 -0.3e1/0.2e1 0.1e1/0.3e1;]/h;
40 S_1=sparse(sparse(1,m)); 50 S_1=sparse(sparse(1,m));
41 S_1(1:4)=S_U; 51 S_1(1:4)=S_U;
42 S_m=sparse(sparse(1,m)); 52 S_m=sparse(sparse(1,m));
43 S_m(m-3:m)=fliplr(-S_U); 53 S_m(m-3:m)=fliplr(-S_U);
44 54