Created
March 5, 2013 14:34
-
-
Save mathsam/5090693 to your computer and use it in GitHub Desktop.
four layer model
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| syms c epsilon1 epsilon2 k N1 N2 gamma | |
| A = [-(-c+0.5*(2*epsilon1+epsilon2+3))*(1+k^2)+1+gamma, -c+0.5*(2*epsilon1+epsilon2+3), 0, 0;... | |
| -c+0.5*(2*epsilon1+epsilon2+1), -(-c+0.5*(2*epsilon1+epsilon2+1))*(k^2+1+N1)+gamma-1+0.5*N1*(1+epsilon1), N1*(-c+0.5*(2*epsilon1+epsilon2+1)),0; ... | |
| 0, N1/epsilon1*(-c+0.5*(epsilon1+epsilon2)), -(-c+0.5*(epsilon1+epsilon2))*(k^2+N1/epsilon1)-N2/epsilon1, N2/epsilon1;... | |
| 0, 0, -c*N2/epsilon2, c*(k^2+N2/epsilon2)+gamma-0.5*N2/epsilon2*(epsilon1+epsilon2)]; | |
| det_A = vpa(det(A)); | |
| %% | |
| clear | |
| wavenumberSet = 0.01:0.02:2; | |
| criticalitySet = 0.01:0.02:4; | |
| c1 = zeros(length(criticalitySet),length(wavenumberSet)); | |
| c2 = zeros(length(criticalitySet),length(wavenumberSet)); | |
| c3 = zeros(length(criticalitySet),length(wavenumberSet)); | |
| c4 = zeros(length(criticalitySet),length(wavenumberSet)); | |
| N1 = 0.2; | |
| N2 = 0.2; | |
| epsilon1 = 5/(2.5+1); | |
| epsilon2 = 5/2; | |
| for i_wavenumber = 1:length(wavenumberSet) | |
| for i_criticality = 1:length(criticalitySet) | |
| k = wavenumberSet(i_wavenumber); | |
| gamma = criticalitySet(i_criticality); | |
| poly_A(1) = (0.0625*(32.*epsilon1*epsilon2*k^6 + 16.*epsilon1*epsilon2*k^8 + 32.*epsilon2*k^4*N1 + 16.*epsilon1*epsilon2*k^4*N1 + 16.*epsilon2*k^6*N1 + 16.*epsilon1*epsilon2*k^6*N1 + 32.*epsilon1*k^4*N2 + 16.*epsilon1*k^6*N2 + 32.*k^2*N1*N2 + 16.*epsilon1*k^2*N1*N2 + 16.*k^4*N1*N2 + 16.*epsilon1*k^4*N1*N2))/(epsilon1*epsilon2); | |
| poly_A(2) = (0.0625*(64.*epsilon1*epsilon2*gamma*k^4 - 64.*epsilon1*epsilon2*k^6 - 80.*epsilon1^2*epsilon2*k^6 - 48.*epsilon1*epsilon2^2*k^6 + 48.*epsilon1*epsilon2*gamma*k^6 - 32.*epsilon1*epsilon2*k^8 - 40.*epsilon1^2*epsilon2*k^8 - 24.*epsilon1*epsilon2^2*k^8 + 64.*epsilon2*gamma*k^2*N1 + 16.*epsilon1*epsilon2*gamma*k^2*N1 - 64.*epsilon2*k^4*N1 - 88.*epsilon1*epsilon2*k^4*N1 - 32.*epsilon1^2*epsilon2*k^4*N1 - 48.*epsilon2^2*k^4*N1 - 24.*epsilon1*epsilon2^2*k^4*N1 + 48.*epsilon2*gamma*k^4*N1 + 32.*epsilon1*epsilon2*gamma*k^4*N1 - 32.*epsilon2*k^6*N1 - 64.*epsilon1*epsilon2*k^6*N1 - 32.*epsilon1^2*epsilon2*k^6*N1 - 24.*epsilon2^2*k^6*N1 - 24.*epsilon1*epsilon2^2*k^6*N1 + 8.*epsilon2*k^2*N1^2 + 8.*epsilon1*epsilon2*k^2*N1^2 + 8.*epsilon2*k^4*N1^2 + 8.*epsilon1*epsilon2*k^4*N1^2 + 32.*epsilon1*gamma*k^2*N2 - 64.*epsilon1*k^4*N2 - 96.*epsilon1^2*k^4*N2 - 32.*epsilon2*k^4*N2 - 64.*epsilon1*epsilon2*k^4*N2 + 32.*epsilon1*gamma*k^4*N2 - 32.*epsilon1*k^6*N2 - 48.*epsilon1^2*k^6*N2 - 16.*epsilon2*k^6*N2 - 32.*epsilon1*epsilon2*k^6*N2 + 32.*gamma*N1*N2 - 64.*k^2*N1*N2 - 104.*epsilon1*k^2*N1*N2 - 40.*epsilon1^2*k^2*N1*N2 - 80.*epsilon2*k^2*N1*N2 - 32.*epsilon1*epsilon2*k^2*N1*N2 + 32.*gamma*k^2*N1*N2 + 16.*epsilon1*gamma*k^2*N1*N2 - 32.*k^4*N1*N2 - 72.*epsilon1*k^4*N1*N2 - 40.*epsilon1^2*k^4*N1*N2 - 48.*epsilon2*k^4*N1*N2 - 32.*epsilon1*epsilon2*k^4*N1*N2 + 8.*N1^2*N2 + 8.*epsilon1*N1^2*N2 + 8.*k^2*N1^2*N2 + 8.*epsilon1*k^2*N1^2*N2))/(epsilon1*epsilon2); | |
| poly_A(3) = (0.0625*(32.*epsilon1*epsilon2*gamma^2*k^2 - 96.*epsilon1*epsilon2*gamma*k^4 - 128.*epsilon1^2*epsilon2*gamma*k^4 - 80.*epsilon1*epsilon2^2*gamma*k^4 + 48.*epsilon1*epsilon2*gamma^2*k^4 + 40.*epsilon1*epsilon2*k^6 + 96.*epsilon1^2*epsilon2*k^6 + 64.*epsilon1^3*epsilon2*k^6 + 64.*epsilon1*epsilon2^2*k^6 + 80.*epsilon1^2*epsilon2^2*k^6 + 24.*epsilon1*epsilon2^3*k^6 - 64.*epsilon1*epsilon2*gamma*k^6 - 88.*epsilon1^2*epsilon2*gamma*k^6 - 56.*epsilon1*epsilon2^2*gamma*k^6 + 12.*epsilon1*epsilon2*k^8 + 48.*epsilon1^2*epsilon2*k^8 + 32.*epsilon1^3*epsilon2*k^8 + 32.*epsilon1*epsilon2^2*k^8 + 40.*epsilon1^2*epsilon2^2*k^8 + 12.*epsilon1*epsilon2^3*k^8 + 32.*epsilon2*gamma^2*N1 - 96.*epsilon2*gamma*k^2*N1 - 136.*epsilon1*epsilon2*gamma*k^2*N1 - 32.*epsilon1^2*epsilon2*gamma*k^2*N1 - 80.*epsilon2^2*gamma*k^2*N1 - 24.*epsilon1*epsilon2^2*gamma*k^2*N1 + 48.*epsilon2*gamma^2*k^2*N1 + 16.*epsilon1*epsilon2*gamma^2*k^2*N1 + 40.*epsilon2*k^4*N1 + 96.*epsilon1*epsilon2*k^4*N1 + 72.*epsilon1^2*epsilon2*k^4*N1 + 20.*epsilon1^3*epsilon2*k^4*N1 + 64.*epsilon2^2*k^4*N1 + 88.*epsilon1*epsilon2^2*k^4*N1 + 32.*epsilon1^2*epsilon2^2*k^4*N1 + 24.*epsilon2^3*k^4*N1 + 12.*epsilon1*epsilon2^3*k^4*N1 - 64.*epsilon2*gamma*k^4*N1 - 112.*epsilon1*epsilon2*gamma*k^4*N1 - 48.*epsilon1^2*epsilon2*gamma*k^4*N1 - 56.*epsilon2^2*gamma*k^4*N1 - 40.*epsilon1*epsilon2^2*gamma*k^4*N1 + 12.*epsilon2*k^6*N1 + 48.*epsilon1*epsilon2*k^6*N1 + 56.*epsilon1^2*epsilon2*k^6*N1 + 20.*epsilon1^3*epsilon2*k^6*N1 + 32.*epsilon2^2*k^6*N1 + 64.*epsilon1*epsilon2^2*k^6*N1 + 32.*epsilon1^2*epsilon2^2*k^6*N1 + 12.*epsilon2^3*k^6*N1 + 12.*epsilon1*epsilon2^3*k^6*N1 + 8.*epsilon2*gamma*N1^2 + 8.*epsilon1*epsilon2*gamma*N1^2 - 4.*epsilon2*k^2*N1^2 - 16.*epsilon1*epsilon2*k^2*N1^2 - 12.*epsilon1^2*epsilon2*k^2*N1^2 - 8.*epsilon2^2*k^2*N1^2 - 8.*epsilon1*epsilon2^2*k^2*N1^2 + 16.*epsilon2*gamma*k^2*N1^2 + 16.*epsilon1*epsilon2*gamma*k^2*N1^2 - 12.*epsilon2*k^4*N1^2 - 24.*epsilon1*epsilon2*k^4*N1^2 - 12.*epsilon1^2*epsilon2*k^4*N1^2 - 8.*epsilon2^2*k^4*N1^2 - 8.*epsilon1*epsilon2^2*k^4*N1^2 - 32.*epsilon1*gamma*k^2*N2 - 64.*epsilon1^2*gamma*k^2*N2 - 64.*epsilon2*gamma*k^2*N2 - 48.*epsilon1*epsilon2*gamma*k^2*N2 + 16.*epsilon1*gamma^2*k^2*N2 + 40.*epsilon1*k^4*N2 + 128.*epsilon1^2*k^4*N2 + 104.*epsilon1^3*k^4*N2 + 64.*epsilon2*k^4*N2 + 160.*epsilon1*epsilon2*k^4*N2 + 144.*epsilon1^2*epsilon2*k^4*N2 + 32.*epsilon2^2*k^4*N2 + 48.*epsilon1*epsilon2^2*k^4*N2 - 32.*epsilon1*gamma*k^4*N2 - 64.*epsilon1^2*gamma*k^4*N2 - 48.*epsilon2*gamma*k^4*N2 - 48.*epsilon1*epsilon2*gamma*k^4*N2 + 12.*epsilon1*k^6*N2 + 64.*epsilon1^2*k^6*N2 + 52.*epsilon1^3*k^6*N2 + 32.*epsilon2*k^6*N2 + 80.*epsilon1*epsilon2*k^6*N2 + 72.*epsilon1^2*epsilon2*k^6*N2 + 16.*epsilon2^2*k^6*N2 + 24.*epsilon1*epsilon2^2*k^6*N2 - 32.*gamma*N1*N2 - 64.*epsilon1*gamma*N1*N2 - 64.*epsilon2*gamma*N1*N2 + 16.*gamma^2*N1*N2 + 40.*k^2*N1*N2 + 128.*epsilon1*k^2*N1*N2 + 116.*epsilon1^2*k^2*N1*N2 + 36.*epsilon1^3*k^2*N1*N2 + 104.*epsilon2*k^2*N1*N2 + 180.*epsilon1*epsilon2*k^2*N1*N2 + 60.*epsilon1^2*epsilon2*k^2*N1*N2 + 64.*epsilon2^2*k^2*N1*N2 + 24.*epsilon1*epsilon2^2*k^2*N1*N2 - 32.*gamma*k^2*N1*N2 - 64.*epsilon1*gamma*k^2*N1*N2 - 24.*epsilon1^2*gamma*k^2*N1*N2 - 80.*epsilon2*gamma*k^2*N1*N2 - 24.*epsilon1*epsilon2*gamma*k^2*N1*N2 + 12.*k^4*N1*N2 + 64.*epsilon1*k^4*N1*N2 + 88.*epsilon1^2*k^4*N1*N2 + 36.*epsilon1^3*k^4*N1*N2 + 72.*epsilon2*k^4*N1*N2 + 132.*epsilon1*epsilon2*k^4*N1*N2 + 60.*epsilon1^2*epsilon2*k^4*N1*N2 + 40.*epsilon2^2*k^4*N1*N2 + 24.*epsilon1*epsilon2^2*k^4*N1*N2 - 4.*N1^2*N2 - 20.*epsilon1*N1^2*N2 - 16.*epsilon1^2*N1^2*N2 - 12.*epsilon2*N1^2*N2 - 12.*epsilon1*epsilon2*N1^2*N2 + 8.*gamma*N1^2*N2 + 8.*epsilon1*gamma*N1^2*N2 - 12.*k^2*N1^2*N2 - 28.*epsilon1*k^2*N1^2*N2 - 16.*epsilon1^2*k^2*N1^2*N2 - 12.*epsilon2*k^2*N1^2*N2 - 12.*epsilon1*epsilon2*k^2*N1^2*N2 + 16.*epsilon1*k^2*N2^2 + 16.*epsilon2*k^2*N2^2 + 8.*epsilon1*k^4*N2^2 + 8.*epsilon2*k^4*N2^2 + 8.*epsilon1*N1*N2^2 + 8.*epsilon2*N1*N2^2 + 8.*epsilon1*k^2*N1*N2^2 + 8.*epsilon2*k^2*N1*N2^2))/(epsilon1*epsilon2); | |
| poly_A(4) = (0.0625*(-32.*epsilon1*epsilon2*gamma^2*k^2 - 48.*epsilon1^2*epsilon2*gamma^2*k^2 - 32.*epsilon1*epsilon2^2*gamma^2*k^2 + 16.*epsilon1*epsilon2*gamma^3*k^2 + 40.*epsilon1*epsilon2*gamma*k^4 + 112.*epsilon1^2*epsilon2*gamma*k^4 + 80.*epsilon1^3*epsilon2*gamma*k^4 + 80.*epsilon1*epsilon2^2*gamma*k^4 + 104.*epsilon1^2*epsilon2^2*gamma*k^4 + 32.*epsilon1*epsilon2^3*gamma*k^4 - 32.*epsilon1*epsilon2*gamma^2*k^4 - 56.*epsilon1^2*epsilon2*gamma^2*k^4 - 40.*epsilon1*epsilon2^2*gamma^2*k^4 - 20.*epsilon1^2*epsilon2*k^6 - 32.*epsilon1^3*epsilon2*k^6 - 16.*epsilon1^4*epsilon2*k^6 - 20.*epsilon1*epsilon2^2*k^6 - 48.*epsilon1^2*epsilon2^2*k^6 - 32.*epsilon1^3*epsilon2^2*k^6 - 16.*epsilon1*epsilon2^3*k^6 - 20.*epsilon1^2*epsilon2^3*k^6 - 4.*epsilon1*epsilon2^4*k^6 + 12.*epsilon1*epsilon2*gamma*k^6 + 64.*epsilon1^2*epsilon2*gamma*k^6 + 48.*epsilon1^3*epsilon2*gamma*k^6 + 48.*epsilon1*epsilon2^2*gamma*k^6 + 64.*epsilon1^2*epsilon2^2*gamma*k^6 + 20.*epsilon1*epsilon2^3*gamma*k^6 - 6.*epsilon1^2*epsilon2*k^8 - 16.*epsilon1^3*epsilon2*k^8 - 8.*epsilon1^4*epsilon2*k^8 - 6.*epsilon1*epsilon2^2*k^8 - 24.*epsilon1^2*epsilon2^2*k^8 - 16.*epsilon1^3*epsilon2^2*k^8 - 8.*epsilon1*epsilon2^3*k^8 - 10.*epsilon1^2*epsilon2^3*k^8 - 2.*epsilon1*epsilon2^4*k^8 - 32.*epsilon2*gamma^2*N1 - 48.*epsilon1*epsilon2*gamma^2*N1 - 32.*epsilon2^2*gamma^2*N1 + 16.*epsilon2*gamma^3*N1 + 40.*epsilon2*gamma*k^2*N1 + 112.*epsilon1*epsilon2*gamma*k^2*N1 + 88.*epsilon1^2*epsilon2*gamma*k^2*N1 + 20.*epsilon1^3*epsilon2*gamma*k^2*N1 + 80.*epsilon2^2*gamma*k^2*N1 + 112.*epsilon1*epsilon2^2*gamma*k^2*N1 + 32.*epsilon1^2*epsilon2^2*gamma*k^2*N1 + 32.*epsilon2^3*gamma*k^2*N1 + 12.*epsilon1*epsilon2^3*gamma*k^2*N1 - 32.*epsilon2*gamma^2*k^2*N1 - 56.*epsilon1*epsilon2*gamma^2*k^2*N1 - 16.*epsilon1^2*epsilon2*gamma^2*k^2*N1 - 40.*epsilon2^2*gamma^2*k^2*N1 - 16.*epsilon1*epsilon2^2*gamma^2*k^2*N1 - 20.*epsilon1*epsilon2*k^4*N1 - 32.*epsilon1^2*epsilon2*k^4*N1 - 18.*epsilon1^3*epsilon2*k^4*N1 - 4.*epsilon1^4*epsilon2*k^4*N1 - 20.*epsilon2^2*k^4*N1 - 48.*epsilon1*epsilon2^2*k^4*N1 - 36.*epsilon1^2*epsilon2^2*k^4*N1 - 10.*epsilon1^3*epsilon2^2*k^4*N1 - 16.*epsilon2^3*k^4*N1 - 22.*epsilon1*epsilon2^3*k^4*N1 - 8.*epsilon1^2*epsilon2^3*k^4*N1 - 4.*epsilon2^4*k^4*N1 - 2.*epsilon1*epsilon2^4*k^4*N1 + 12.*epsilon2*gamma*k^4*N1 + 64.*epsilon1*epsilon2*gamma*k^4*N1 + 72.*epsilon1^2*epsilon2*gamma*k^4*N1 + 24.*epsilon1^3*epsilon2*gamma*k^4*N1 + 48.*epsilon2^2*gamma*k^4*N1 + 88.*epsilon1*epsilon2^2*gamma*k^4*N1 + 40.*epsilon1^2*epsilon2^2*gamma*k^4*N1 + 20.*epsilon2^3*gamma*k^4*N1 + 16.*epsilon1*epsilon2^3*gamma*k^4*N1 - 6.*epsilon1*epsilon2*k^6*N1 - 16.*epsilon1^2*epsilon2*k^6*N1 - 14.*epsilon1^3*epsilon2*k^6*N1 - 4.*epsilon1^4*epsilon2*k^6*N1 - 6.*epsilon2^2*k^6*N1 - 24.*epsilon1*epsilon2^2*k^6*N1 - 28.*epsilon1^2*epsilon2^2*k^6*N1 - 10.*epsilon1^3*epsilon2^2*k^6*N1 - 8.*epsilon2^3*k^6*N1 - 16.*epsilon1*epsilon2^3*k^6*N1 - 8.*epsilon1^2*epsilon2^3*k^6*N1 - 2.*epsilon2^4*k^6*N1 - 2.*epsilon1*epsilon2^4*k^6*N1 - 4.*epsilon2*gamma*N1^2 - 16.*epsilon1*epsilon2*gamma*N1^2 - 12.*epsilon1^2*epsilon2*gamma*N1^2 - 8.*epsilon2^2*gamma*N1^2 - 8.*epsilon1*epsilon2^2*gamma*N1^2 + 8.*epsilon2*gamma^2*N1^2 + 8.*epsilon1*epsilon2*gamma^2*N1^2 + 2.*epsilon1*epsilon2*k^2*N1^2 + 6.*epsilon1^2*epsilon2*k^2*N1^2 + 4.*epsilon1^3*epsilon2*k^2*N1^2 + 2.*epsilon2^2*k^2*N1^2 + 8.*epsilon1*epsilon2^2*k^2*N1^2 + 6.*epsilon1^2*epsilon2^2*k^2*N1^2 + 2.*epsilon2^3*k^2*N1^2 + 2.*epsilon1*epsilon2^3*k^2*N1^2 - 12.*epsilon2*gamma*k^2*N1^2 - 28.*epsilon1*epsilon2*gamma*k^2*N1^2 - 16.*epsilon1^2*epsilon2*gamma*k^2*N1^2 - 12.*epsilon2^2*gamma*k^2*N1^2 - 12.*epsilon1*epsilon2^2*gamma*k^2*N1^2 + 6.*epsilon1*epsilon2*k^4*N1^2 + 10.*epsilon1^2*epsilon2*k^4*N1^2 + 4.*epsilon1^3*epsilon2*k^4*N1^2 + 6.*epsilon2^2*k^4*N1^2 + 12.*epsilon1*epsilon2^2*k^4*N1^2 + 6.*epsilon1^2*epsilon2^2*k^4*N1^2 + 2.*epsilon2^3*k^4*N1^2 + 2.*epsilon1*epsilon2^3*k^4*N1^2 - 32.*epsilon2*gamma^2*N2 + 32.*epsilon1^2*gamma*k^2*N2 + 40.*epsilon1^3*gamma*k^2*N2 + 96.*epsilon2*gamma*k^2*N2 + 128.*epsilon1*epsilon2*gamma*k^2*N2 + 64.*epsilon1^2*epsilon2*gamma*k^2*N2 + 48.*epsilon2^2*gamma*k^2*N2 + 24.*epsilon1*epsilon2^2*gamma*k^2*N2 - 16.*epsilon1^2*gamma^2*k^2*N2 - 48.*epsilon2*gamma^2*k^2*N2 - 16.*epsilon1*epsilon2*gamma^2*k^2*N2 - 40.*epsilon1^2*k^4*N2 - 80.*epsilon1^3*k^4*N2 - 48.*epsilon1^4*k^4*N2 - 40.*epsilon2*k^4*N2 - 104.*epsilon1*epsilon2*k^4*N2 - 160.*epsilon1^2*epsilon2*k^4*N2 - 104.*epsilon1^3*epsilon2*k^4*N2 - 32.*epsilon2^2*k^4*N2 - 80.*epsilon1*epsilon2^2*k^4*N2 - 72.*epsilon1^2*epsilon2^2*k^4*N2 - 8.*epsilon2^3*k^4*N2 - 16.*epsilon1*epsilon2^3*k^4*N2 + 32.*epsilon1^2*gamma*k^4*N2 + 40.*epsilon1^3*gamma*k^4*N2 + 64.*epsilon2*gamma*k^4*N2 + 96.*epsilon1*epsilon2*gamma*k^4*N2 + 64.*epsilon1^2*epsilon2*gamma*k^4*N2 + 32.*epsilon2^2*gamma*k^4*N2 + 24.*epsilon1*epsilon2^2*gamma*k^4*N2 - 12.*epsilon1^2*k^6*N2 - 40.*epsilon1^3*k^6*N2 - 24.*epsilon1^4*k^6*N2 - 12.*epsilon2*k^6*N2 - 44.*epsilon1*epsilon2*k^6*N2 - 80.*epsilon1^2*epsilon2*k^6*N2 - 52.*epsilon1^3*epsilon2*k^6*N2 - 16.*epsilon2^2*k^6*N2 - 40.*epsilon1*epsilon2^2*k^6*N2 - 36.*epsilon1^2*epsilon2^2*k^6*N2 - 4.*epsilon2^3*k^6*N2 - 8.*epsilon1*epsilon2^3*k^6*N2 + 32.*epsilon1*gamma*N1*N2 + 40.*epsilon1^2*gamma*N1*N2 + 40.*epsilon2*gamma*N1*N2 + 88.*epsilon1*epsilon2*gamma*N1*N2 + 40.*epsilon2^2*gamma*N1*N2 - 16.*epsilon1*gamma^2*N1*N2 - 32.*epsilon2*gamma^2*N1*N2 - 40.*epsilon1*k^2*N1*N2 - 80.*epsilon1^2*k^2*N1*N2 - 54.*epsilon1^3*k^2*N1*N2 - 14.*epsilon1^4*k^2*N1*N2 - 40.*epsilon2*k^2*N1*N2 - 132.*epsilon1*epsilon2*k^2*N1*N2 - 124.*epsilon1^2*epsilon2*k^2*N1*N2 - 36.*epsilon1^3*epsilon2*k^2*N1*N2 - 52.*epsilon2^2*k^2*N1*N2 - 90.*epsilon1*epsilon2^2*k^2*N1*N2 - 30.*epsilon1^2*epsilon2^2*k^2*N1*N2 - 20.*epsilon2^3*k^2*N1*N2 - 8.*epsilon1*epsilon2^3*k^2*N1*N2 + 32.*epsilon1*gamma*k^2*N1*N2 + 40.*epsilon1^2*gamma*k^2*N1*N2 + 12.*epsilon1^3*gamma*k^2*N1*N2 + 56.*epsilon2*gamma*k^2*N1*N2 + 96.*epsilon1*epsilon2*gamma*k^2*N1*N2 + 24.*epsilon1^2*epsilon2*gamma*k^2*N1*N2 + 48.*epsilon2^2*gamma*k^2*N1*N2 + 12.*epsilon1*epsilon2^2*gamma*k^2*N1*N2 - 12.*epsilon1*k^4*N1*N2 - 40.*epsilon1^2*k^4*N1*N2 - 42.*epsilon1^3*k^4*N1*N2 - 14.*epsilon1^4*k^4*N1*N2 - 12.*epsilon2*k^4*N1*N2 - 76.*epsilon1*epsilon2*k^4*N1*N2 - 96.*epsilon1^2*epsilon2*k^4*N1*N2 - 36.*epsilon1^3*epsilon2*k^4*N1*N2 - 36.*epsilon2^2*k^4*N1*N2 - 66.*epsilon1*epsilon2^2*k^4*N1*N2 - 30.*epsilon1^2*epsilon2^2*k^4*N1*N2 - 12.*epsilon2^3*k^4*N1*N2 - 8.*epsilon1*epsilon2^3*k^4*N1*N2 + 4.*epsilon1*N1^2*N2 + 14.*epsilon1^2*N1^2*N2 + 10.*epsilon1^3*N1^2*N2 + 4.*epsilon2*N1^2*N2 + 20.*epsilon1*epsilon2*N1^2*N2 + 16.*epsilon1^2*epsilon2*N1^2*N2 + 6.*epsilon2^2*N1^2*N2 + 6.*epsilon1*epsilon2^2*N1^2*N2 - 8.*epsilon1*gamma*N1^2*N2 - 8.*epsilon1^2*gamma*N1^2*N2 - 8.*epsilon2*gamma*N1^2*N2 - 8.*epsilon1*epsilon2*gamma*N1^2*N2 + 12.*epsilon1*k^2*N1^2*N2 + 22.*epsilon1^2*k^2*N1^2*N2 + 10.*epsilon1^3*k^2*N1^2*N2 + 12.*epsilon2*k^2*N1^2*N2 + 28.*epsilon1*epsilon2*k^2*N1^2*N2 + 16.*epsilon1^2*epsilon2*k^2*N1^2*N2 + 6.*epsilon2^2*k^2*N1^2*N2 + 6.*epsilon1*epsilon2^2*k^2*N1^2*N2 + 16.*epsilon1*gamma*N2^2 + 16.*epsilon2*gamma*N2^2 - 32.*epsilon1*k^2*N2^2 - 32.*epsilon1^2*k^2*N2^2 - 32.*epsilon2*k^2*N2^2 - 48.*epsilon1*epsilon2*k^2*N2^2 - 16.*epsilon2^2*k^2*N2^2 + 16.*epsilon1*gamma*k^2*N2^2 + 16.*epsilon2*gamma*k^2*N2^2 - 16.*epsilon1*k^4*N2^2 - 16.*epsilon1^2*k^4*N2^2 - 16.*epsilon2*k^4*N2^2 - 24.*epsilon1*epsilon2*k^4*N2^2 - 8.*epsilon2^2*k^4*N2^2 - 4.*epsilon1*N1*N2^2 - 12.*epsilon1^2*N1*N2^2 - 4.*epsilon2*N1*N2^2 - 20.*epsilon1*epsilon2*N1*N2^2 - 8.*epsilon2^2*N1*N2^2 + 8.*epsilon1*gamma*N1*N2^2 + 8.*epsilon2*gamma*N1*N2^2 - 12.*epsilon1*k^2*N1*N2^2 - 12.*epsilon1^2*k^2*N1*N2^2 - 12.*epsilon2*k^2*N1*N2^2 - 20.*epsilon1*epsilon2*k^2*N1*N2^2 - 8.*epsilon2^2*k^2*N1*N2^2))/(epsilon1*epsilon2); | |
| poly_A(5) = (0.0625*(16.*epsilon1^2*epsilon2*gamma^2*k^2 + 16.*epsilon1^3*epsilon2*gamma^2*k^2 + 16.*epsilon1*epsilon2^2*gamma^2*k^2 + 24.*epsilon1^2*epsilon2^2*gamma^2*k^2 + 8.*epsilon1*epsilon2^3*gamma^2*k^2 - 8.*epsilon1^2*epsilon2*gamma^3*k^2 - 8.*epsilon1*epsilon2^2*gamma^3*k^2 - 20.*epsilon1^2*epsilon2*gamma*k^4 - 32.*epsilon1^3*epsilon2*gamma*k^4 - 16.*epsilon1^4*epsilon2*gamma*k^4 - 20.*epsilon1*epsilon2^2*gamma*k^4 - 48.*epsilon1^2*epsilon2^2*gamma*k^4 - 32.*epsilon1^3*epsilon2^2*gamma*k^4 - 16.*epsilon1*epsilon2^3*gamma*k^4 - 20.*epsilon1^2*epsilon2^3*gamma*k^4 - 4.*epsilon1*epsilon2^4*gamma*k^4 + 16.*epsilon1^2*epsilon2*gamma^2*k^4 + 16.*epsilon1^3*epsilon2*gamma^2*k^4 + 16.*epsilon1*epsilon2^2*gamma^2*k^4 + 24.*epsilon1^2*epsilon2^2*gamma^2*k^4 + 8.*epsilon1*epsilon2^3*gamma^2*k^4 - 6.*epsilon1^2*epsilon2*gamma*k^6 - 16.*epsilon1^3*epsilon2*gamma*k^6 - 8.*epsilon1^4*epsilon2*gamma*k^6 - 6.*epsilon1*epsilon2^2*gamma*k^6 - 24.*epsilon1^2*epsilon2^2*gamma*k^6 - 16.*epsilon1^3*epsilon2^2*gamma*k^6 - 8.*epsilon1*epsilon2^3*gamma*k^6 - 10.*epsilon1^2*epsilon2^3*gamma*k^6 - 2.*epsilon1*epsilon2^4*gamma*k^6 + 16.*epsilon1*epsilon2*gamma^2*N1 + 16.*epsilon1^2*epsilon2*gamma^2*N1 + 16.*epsilon2^2*gamma^2*N1 + 24.*epsilon1*epsilon2^2*gamma^2*N1 + 8.*epsilon2^3*gamma^2*N1 - 8.*epsilon1*epsilon2*gamma^3*N1 - 8.*epsilon2^2*gamma^3*N1 - 20.*epsilon1*epsilon2*gamma*k^2*N1 - 32.*epsilon1^2*epsilon2*gamma*k^2*N1 - 18.*epsilon1^3*epsilon2*gamma*k^2*N1 - 4.*epsilon1^4*epsilon2*gamma*k^2*N1 - 20.*epsilon2^2*gamma*k^2*N1 - 48.*epsilon1*epsilon2^2*gamma*k^2*N1 - 36.*epsilon1^2*epsilon2^2*gamma*k^2*N1 - 10.*epsilon1^3*epsilon2^2*gamma*k^2*N1 - 16.*epsilon2^3*gamma*k^2*N1 - 22.*epsilon1*epsilon2^3*gamma*k^2*N1 - 8.*epsilon1^2*epsilon2^3*gamma*k^2*N1 - 4.*epsilon2^4*gamma*k^2*N1 - 2.*epsilon1*epsilon2^4*gamma*k^2*N1 + 16.*epsilon1*epsilon2*gamma^2*k^2*N1 + 16.*epsilon1^2*epsilon2*gamma^2*k^2*N1 + 4.*epsilon1^3*epsilon2*gamma^2*k^2*N1 + 16.*epsilon2^2*gamma^2*k^2*N1 + 24.*epsilon1*epsilon2^2*gamma^2*k^2*N1 + 8.*epsilon1^2*epsilon2^2*gamma^2*k^2*N1 + 8.*epsilon2^3*gamma^2*k^2*N1 + 4.*epsilon1*epsilon2^3*gamma^2*k^2*N1 - 6.*epsilon1*epsilon2*gamma*k^4*N1 - 16.*epsilon1^2*epsilon2*gamma*k^4*N1 - 14.*epsilon1^3*epsilon2*gamma*k^4*N1 - 4.*epsilon1^4*epsilon2*gamma*k^4*N1 - 6.*epsilon2^2*gamma*k^4*N1 - 24.*epsilon1*epsilon2^2*gamma*k^4*N1 - 28.*epsilon1^2*epsilon2^2*gamma*k^4*N1 - 10.*epsilon1^3*epsilon2^2*gamma*k^4*N1 - 8.*epsilon2^3*gamma*k^4*N1 - 16.*epsilon1*epsilon2^3*gamma*k^4*N1 - 8.*epsilon1^2*epsilon2^3*gamma*k^4*N1 - 2.*epsilon2^4*gamma*k^4*N1 - 2.*epsilon1*epsilon2^4*gamma*k^4*N1 + 2.*epsilon1*epsilon2*gamma*N1^2 + 6.*epsilon1^2*epsilon2*gamma*N1^2 + 4.*epsilon1^3*epsilon2*gamma*N1^2 + 2.*epsilon2^2*gamma*N1^2 + 8.*epsilon1*epsilon2^2*gamma*N1^2 + 6.*epsilon1^2*epsilon2^2*gamma*N1^2 + 2.*epsilon2^3*gamma*N1^2 + 2.*epsilon1*epsilon2^3*gamma*N1^2 - 4.*epsilon1*epsilon2*gamma^2*N1^2 - 4.*epsilon1^2*epsilon2*gamma^2*N1^2 - 4.*epsilon2^2*gamma^2*N1^2 - 4.*epsilon1*epsilon2^2*gamma^2*N1^2 + 6.*epsilon1*epsilon2*gamma*k^2*N1^2 + 10.*epsilon1^2*epsilon2*gamma*k^2*N1^2 + 4.*epsilon1^3*epsilon2*gamma*k^2*N1^2 + 6.*epsilon2^2*gamma*k^2*N1^2 + 12.*epsilon1*epsilon2^2*gamma*k^2*N1^2 + 6.*epsilon1^2*epsilon2^2*gamma*k^2*N1^2 + 2.*epsilon2^3*gamma*k^2*N1^2 + 2.*epsilon1*epsilon2^3*gamma*k^2*N1^2 + 32.*epsilon2*gamma^2*N2 + 32.*epsilon1*epsilon2*gamma^2*N2 + 16.*epsilon2^2*gamma^2*N2 - 16.*epsilon2*gamma^3*N2 - 8.*epsilon1^3*gamma*k^2*N2 - 8.*epsilon1^4*gamma*k^2*N2 - 40.*epsilon2*gamma*k^2*N2 - 64.*epsilon1*epsilon2*gamma*k^2*N2 - 48.*epsilon1^2*epsilon2*gamma*k^2*N2 - 20.*epsilon1^3*epsilon2*gamma*k^2*N2 - 32.*epsilon2^2*gamma*k^2*N2 - 40.*epsilon1*epsilon2^2*gamma*k^2*N2 - 16.*epsilon1^2*epsilon2^2*gamma*k^2*N2 - 8.*epsilon2^3*gamma*k^2*N2 - 4.*epsilon1*epsilon2^3*gamma*k^2*N2 + 4.*epsilon1^3*gamma^2*k^2*N2 + 32.*epsilon2*gamma^2*k^2*N2 + 32.*epsilon1*epsilon2*gamma^2*k^2*N2 + 8.*epsilon1^2*epsilon2*gamma^2*k^2*N2 + 16.*epsilon2^2*gamma^2*k^2*N2 + 4.*epsilon1*epsilon2^2*gamma^2*k^2*N2 + 10.*epsilon1^3*k^4*N2 + 16.*epsilon1^4*k^4*N2 + 8.*epsilon1^5*k^4*N2 + 20.*epsilon1^2*epsilon2*k^4*N2 + 40.*epsilon1^3*epsilon2*k^4*N2 + 24.*epsilon1^4*epsilon2*k^4*N2 + 10.*epsilon1*epsilon2^2*k^4*N2 + 32.*epsilon1^2*epsilon2^2*k^4*N2 + 26.*epsilon1^3*epsilon2^2*k^4*N2 + 8.*epsilon1*epsilon2^3*k^4*N2 + 12.*epsilon1^2*epsilon2^3*k^4*N2 + 2.*epsilon1*epsilon2^4*k^4*N2 - 8.*epsilon1^3*gamma*k^4*N2 - 8.*epsilon1^4*gamma*k^4*N2 - 12.*epsilon2*gamma*k^4*N2 - 32.*epsilon1*epsilon2*gamma*k^4*N2 - 32.*epsilon1^2*epsilon2*gamma*k^4*N2 - 20.*epsilon1^3*epsilon2*gamma*k^4*N2 - 16.*epsilon2^2*gamma*k^4*N2 - 24.*epsilon1*epsilon2^2*gamma*k^4*N2 - 16.*epsilon1^2*epsilon2^2*gamma*k^4*N2 - 4.*epsilon2^3*gamma*k^4*N2 - 4.*epsilon1*epsilon2^3*gamma*k^4*N2 + 3.*epsilon1^3*k^6*N2 + 8.*epsilon1^4*k^6*N2 + 4.*epsilon1^5*k^6*N2 + 6.*epsilon1^2*epsilon2*k^6*N2 + 20.*epsilon1^3*epsilon2*k^6*N2 + 12.*epsilon1^4*epsilon2*k^6*N2 + 3.*epsilon1*epsilon2^2*k^6*N2 + 16.*epsilon1^2*epsilon2^2*k^6*N2 + 13.*epsilon1^3*epsilon2^2*k^6*N2 + 4.*epsilon1*epsilon2^3*k^6*N2 + 6.*epsilon1^2*epsilon2^3*k^6*N2 + epsilon1*epsilon2^4*k^6*N2 - 8.*epsilon1^2*gamma*N1*N2 - 8.*epsilon1^3*gamma*N1*N2 - 20.*epsilon1*epsilon2*gamma*N1*N2 - 28.*epsilon1^2*epsilon2*gamma*N1*N2 - 12.*epsilon2^2*gamma*N1*N2 - 28.*epsilon1*epsilon2^2*gamma*N1*N2 - 8.*epsilon2^3*gamma*N1*N2 + 4.*epsilon1^2*gamma^2*N1*N2 + 16.*epsilon1*epsilon2*gamma^2*N1*N2 + 12.*epsilon2^2*gamma^2*N1*N2 + 10.*epsilon1^2*k^2*N1*N2 + 16.*epsilon1^3*k^2*N1*N2 + 9.*epsilon1^4*k^2*N1*N2 + 2.*epsilon1^5*k^2*N1*N2 + 20.*epsilon1*epsilon2*k^2*N1*N2 + 40.*epsilon1^2*epsilon2*k^2*N1*N2 + 27.*epsilon1^3*epsilon2*k^2*N1*N2 + 7.*epsilon1^4*epsilon2*k^2*N1*N2 + 10.*epsilon2^2*k^2*N1*N2 + 32.*epsilon1*epsilon2^2*k^2*N1*N2 + 29.*epsilon1^2*epsilon2^2*k^2*N1*N2 + 9.*epsilon1^3*epsilon2^2*k^2*N1*N2 + 8.*epsilon2^3*k^2*N1*N2 + 13.*epsilon1*epsilon2^3*k^2*N1*N2 + 5.*epsilon1^2*epsilon2^3*k^2*N1*N2 + 2.*epsilon2^4*k^2*N1*N2 + epsilon1*epsilon2^4*k^2*N1*N2 - 8.*epsilon1^2*gamma*k^2*N1*N2 - 8.*epsilon1^3*gamma*k^2*N1*N2 - 2.*epsilon1^4*gamma*k^2*N1*N2 - 28.*epsilon1*epsilon2*gamma*k^2*N1*N2 - 28.*epsilon1^2*epsilon2*gamma*k^2*N1*N2 - 6.*epsilon1^3*epsilon2*gamma*k^2*N1*N2 - 20.*epsilon2^2*gamma*k^2*N1*N2 - 28.*epsilon1*epsilon2^2*gamma*k^2*N1*N2 - 6.*epsilon1^2*epsilon2^2*gamma*k^2*N1*N2 - 8.*epsilon2^3*gamma*k^2*N1*N2 - 2.*epsilon1*epsilon2^3*gamma*k^2*N1*N2 + 3.*epsilon1^2*k^4*N1*N2 + 8.*epsilon1^3*k^4*N1*N2 + 7.*epsilon1^4*k^4*N1*N2 + 2.*epsilon1^5*k^4*N1*N2 + 6.*epsilon1*epsilon2*k^4*N1*N2 + 20.*epsilon1^2*epsilon2*k^4*N1*N2 + 21.*epsilon1^3*epsilon2*k^4*N1*N2 + 7.*epsilon1^4*epsilon2*k^4*N1*N2 + 3.*epsilon2^2*k^4*N1*N2 + 16.*epsilon1*epsilon2^2*k^4*N1*N2 + 22.*epsilon1^2*epsilon2^2*k^4*N1*N2 + 9.*epsilon1^3*epsilon2^2*k^4*N1*N2 + 4.*epsilon2^3*k^4*N1*N2 + 9.*epsilon1*epsilon2^3*k^4*N1*N2 + 5.*epsilon1^2*epsilon2^3*k^4*N1*N2 + epsilon2^4*k^4*N1*N2 + epsilon1*epsilon2^4*k^4*N1*N2 - 1.*epsilon1^2*N1^2*N2 - 3.*epsilon1^3*N1^2*N2 - 2.*epsilon1^4*N1^2*N2 - 2.*epsilon1*epsilon2*N1^2*N2 - 7.*epsilon1^2*epsilon2*N1^2*N2 - 5.*epsilon1^3*epsilon2*N1^2*N2 - 1.*epsilon2^2*N1^2*N2 - 5.*epsilon1*epsilon2^2*N1^2*N2 - 4.*epsilon1^2*epsilon2^2*N1^2*N2 - 1.*epsilon2^3*N1^2*N2 - 1.*epsilon1*epsilon2^3*N1^2*N2 + 2.*epsilon1^2*gamma*N1^2*N2 + 2.*epsilon1^3*gamma*N1^2*N2 + 4.*epsilon1*epsilon2*gamma*N1^2*N2 + 4.*epsilon1^2*epsilon2*gamma*N1^2*N2 + 2.*epsilon2^2*gamma*N1^2*N2 + 2.*epsilon1*epsilon2^2*gamma*N1^2*N2 - 3.*epsilon1^2*k^2*N1^2*N2 - 5.*epsilon1^3*k^2*N1^2*N2 - 2.*epsilon1^4*k^2*N1^2*N2 - 6.*epsilon1*epsilon2*k^2*N1^2*N2 - 11.*epsilon1^2*epsilon2*k^2*N1^2*N2 - 5.*epsilon1^3*epsilon2*k^2*N1^2*N2 - 3.*epsilon2^2*k^2*N1^2*N2 - 7.*epsilon1*epsilon2^2*k^2*N1^2*N2 - 4.*epsilon1^2*epsilon2^2*k^2*N1^2*N2 - 1.*epsilon2^3*k^2*N1^2*N2 - 1.*epsilon1*epsilon2^3*k^2*N1^2*N2 - 16.*epsilon1*gamma*N2^2 - 16.*epsilon1^2*gamma*N2^2 - 16.*epsilon2*gamma*N2^2 - 24.*epsilon1*epsilon2*gamma*N2^2 - 8.*epsilon2^2*gamma*N2^2 + 8.*epsilon1*gamma^2*N2^2 + 8.*epsilon2*gamma^2*N2^2 + 20.*epsilon1*k^2*N2^2 + 32.*epsilon1^2*k^2*N2^2 + 16.*epsilon1^3*k^2*N2^2 + 20.*epsilon2*k^2*N2^2 + 48.*epsilon1*epsilon2*k^2*N2^2 + 32.*epsilon1^2*epsilon2*k^2*N2^2 + 16.*epsilon2^2*k^2*N2^2 + 20.*epsilon1*epsilon2^2*k^2*N2^2 + 4.*epsilon2^3*k^2*N2^2 - 16.*epsilon1*gamma*k^2*N2^2 - 16.*epsilon1^2*gamma*k^2*N2^2 - 16.*epsilon2*gamma*k^2*N2^2 - 24.*epsilon1*epsilon2*gamma*k^2*N2^2 - 8.*epsilon2^2*gamma*k^2*N2^2 + 6.*epsilon1*k^4*N2^2 + 16.*epsilon1^2*k^4*N2^2 + 8.*epsilon1^3*k^4*N2^2 + 6.*epsilon2*k^4*N2^2 + 24.*epsilon1*epsilon2*k^4*N2^2 + 16.*epsilon1^2*epsilon2*k^4*N2^2 + 8.*epsilon2^2*k^4*N2^2 + 10.*epsilon1*epsilon2^2*k^4*N2^2 + 2.*epsilon2^3*k^4*N2^2 + 2.*epsilon1^2*N1*N2^2 + 4.*epsilon1^3*N1*N2^2 + 4.*epsilon1*epsilon2*N1*N2^2 + 10.*epsilon1^2*epsilon2*N1*N2^2 + 2.*epsilon2^2*N1*N2^2 + 8.*epsilon1*epsilon2^2*N1*N2^2 + 2.*epsilon2^3*N1*N2^2 - 4.*epsilon1^2*gamma*N1*N2^2 - 8.*epsilon1*epsilon2*gamma*N1*N2^2 - 4.*epsilon2^2*gamma*N1*N2^2 + 6.*epsilon1^2*k^2*N1*N2^2 + 4.*epsilon1^3*k^2*N1*N2^2 + 12.*epsilon1*epsilon2*k^2*N1*N2^2 + 10.*epsilon1^2*epsilon2*k^2*N1*N2^2 + 6.*epsilon2^2*k^2*N1*N2^2 + 8.*epsilon1*epsilon2^2*k^2*N1*N2^2 + 2.*epsilon2^3*k^2*N1*N2^2))/(epsilon1*epsilon2); | |
| waveSpeed = roots(poly_A); | |
| c1(i_criticality,i_wavenumber) = waveSpeed(1); | |
| c2(i_criticality,i_wavenumber) = waveSpeed(2); | |
| c3(i_criticality,i_wavenumber) = waveSpeed(3); | |
| c4(i_criticality,i_wavenumber) = waveSpeed(4); | |
| end | |
| end | |
| %% plot | |
| wavenumber2d = repmat(wavenumberSet(:)',length(criticalitySet),1); | |
| tmp1 = imag(c1.*wavenumber2d); | |
| tmp2 = imag(c2.*wavenumber2d); | |
| tmp3 = imag(c3.*wavenumber2d); | |
| tmp4 = imag(c4.*wavenumber2d); | |
| growthRate1 = zeros(size(tmp1)); | |
| growthRate2 = zeros(size(tmp1)); | |
| tmpindex = tmp1(:)>0; | |
| growthRate1(tmpindex) = tmp1(tmpindex); | |
| tmpindex = tmp2(:)>0; | |
| growthRate1(tmpindex) = tmp2(tmpindex); | |
| tmpindex = tmp3(:)>0; | |
| growthRate1(tmpindex) = tmp3(tmpindex); | |
| tmpindex = tmp4(:)>0; | |
| growthRate1(tmpindex) = tmp4(tmpindex); | |
| tmpindex = tmp1(:)<0; | |
| growthRate2(tmpindex) = tmp1(tmpindex); | |
| tmpindex = tmp2(:)<0; | |
| growthRate2(tmpindex) = tmp2(tmpindex); | |
| tmpindex = tmp3(:)<0; | |
| growthRate2(tmpindex) = tmp3(tmpindex); | |
| tmpindex = tmp4(:)<0; | |
| growthRate2(tmpindex) = tmp4(tmpindex); | |
| figure | |
| contourf(wavenumberSet,criticalitySet,growthRate1,30) | |
| colorbar | |
| set(gca,'YDir','reverse') | |
| xlabel('wavenumber','fontsize',16) | |
| ylabel('1/Sc','fontsize',16) | |
| set(gca,'fontsize',16) | |
| xlim([0 2]) | |
| %{ | |
| figure | |
| contourf(wavenumberSet,criticalitySet,growthRate2,30) | |
| colorbar | |
| set(gca,'YDir','reverse') | |
| xlabel('wavenumber','fontsize',16) | |
| ylabel('1/Sc','fontsize',16) | |
| set(gca,'fontsize',16) | |
| xlim([0 2]) | |
| %ylim([0.2 0.5]) | |
| %} | |
| %% calculate the eigen vectors correspond to phase speed (growthRate>0) | |
| wavenumber2d = repmat(wavenumberSet(:)',length(criticalitySet),1); | |
| tmp1 = c1.*wavenumber2d; | |
| tmp2 = c2.*wavenumber2d; | |
| tmp3 = c3.*wavenumber2d; | |
| tmp4 = c4.*wavenumber2d; | |
| mod1 = zeros(size(tmp1)); | |
| tmpindex = imag(tmp1(:))>0; | |
| mod1(tmpindex) = tmp1(tmpindex); | |
| tmpindex = imag(tmp2(:))>0; | |
| mod1(tmpindex) = tmp2(tmpindex); | |
| tmpindex = imag(tmp3(:))>0; | |
| mod1(tmpindex) = tmp3(tmpindex); | |
| tmpindex = imag(tmp4(:))>0; | |
| mod1(tmpindex) = tmp4(tmpindex); | |
| phi = zeros(size(c1,1),size(c1,2),3); | |
| for i_wavenumber = 1:length(wavenumberSet) | |
| for i_criticality = 1:length(criticalitySet) | |
| c = mod1(i_criticality,i_wavenumber); | |
| k = wavenumberSet(i_wavenumber); | |
| gamma = criticalitySet(i_criticality); | |
| A = [-(-c+0.5*(2*epsilon1+epsilon2+3))*(1+k^2)+1+gamma, -c+0.5*(2*epsilon1+epsilon2+3), 0, 0;... | |
| -c+0.5*(2*epsilon1+epsilon2+1), -(-c+0.5*(2*epsilon1+epsilon2+1))*(k^2+1+N1)+gamma-1+0.5*N1*(1+epsilon1), N1*(-c+0.5*(2*epsilon1+epsilon2+1)),0; ... | |
| 0, N1/epsilon1*(-c+0.5*(epsilon1+epsilon2)), -(-c+0.5*(epsilon1+epsilon2))*(k^2+N1/epsilon1)-N2/epsilon1, N2/epsilon1;... | |
| 0, 0, -c*N2/epsilon2, c*(k^2+N2/epsilon2)+gamma-0.5*N2/epsilon2*(epsilon1+epsilon2)]; | |
| Areduced = [A(1,1),A(1,2),A(1,3);... | |
| A(2,1),A(2,2),A(2,3);... | |
| A(3,1),A(3,2),A(3,3)]; | |
| b = -[A(1,4);A(2,4);A(3,4)]; | |
| phi(i_criticality,i_wavenumber,:) = Areduced\b; | |
| end | |
| end | |
| %% | |
| mask_growthRate = zeros(size(growthRate1)); | |
| tmpindex = growthRate1>0; | |
| mask_growthRate(tmpindex) = 1; | |
| figure | |
| subplot(3,1,1) | |
| imagesc(wavenumberSet,criticalitySet,real((phi(:,:,1))).*mask_growthRate) | |
| caxis([-1 1]) | |
| hold on | |
| [C1,H1]=contour(wavenumberSet,criticalitySet,real((phi(:,:,1))).*mask_growthRate,[-1 -1],'w-'); | |
| set(gca,'YDir','reverse') | |
| subplot(3,1,2) | |
| imagesc(wavenumberSet,criticalitySet,real((phi(:,:,2))).*mask_growthRate) | |
| caxis([-1 1]) | |
| set(gca,'YDir','reverse') | |
| subplot(3,1,3) | |
| imagesc(wavenumberSet,criticalitySet,real((phi(:,:,2))).*mask_growthRate) | |
| caxis([-1 1]) | |
| set(gca,'YDir','reverse') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment