Kỹ thuật mô phỏng quang - Quang phổ và vật lý Plasma

58 trang phuongnguyen 4760

Download

Bạn đang xem 20 trang mẫu của tài liệu "Kỹ thuật mô phỏng quang - Quang phổ và vật lý Plasma", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

ky_thuat_mo_phong_quang_quang_pho_va_vat_ly_plasma.doc

Nội dung text: Kỹ thuật mô phỏng quang - Quang phổ và vật lý Plasma

SWEET NOVEMBER 1 VERSION 2009  Kỹ thuật mô phỏng quang - quang phổ và vật lý Plasma VIC FAN
SWEET NOVEMBER 2 VERSION 2009 (Bài viết này dùng cho phiên bản Matlab 6.5. Bạn đọc cũng nên sử dụng phiên bản này vì kích thước gọn nhẹ, đồng thời vẫn đảm bảo các tính năng cần thiết cho bài mô phỏng) A. NHẬP MÔN MATLAB 1. Cửa sổ làm việc  Sau khi cài đặt phần mềm thành công, lần đầu tiên khởi động Matlab, giao diện chương trình sẽ xuất hiện gồm nhiều khung làm việc (Workspace, Command History, Direct History, Command Window, ). Bạn hãy tắt tất cả chúng, chỉ giữ lại Command Window. Workspace Command Window Bạn sẽ viết phần lập trình của mình (thuật ngữ gọi là “code”) trong cửa sổ M-file và chạy chương trình trong Command Window (bằng Command History cách nhấn menu Debug > Save and Run (hay Run)) Cửa sổ M-file  Vào menu File > New > M-file, lúc này sẽ có một cửa sổ mới xuất hiện. Code (Phần chương trình viết) 2. Một vài toán tử so sánh VIC FAN
SWEET NOVEMBER 3 VERSION 2009 Toán tử quan hệ Ý nghĩa > Lớn hơn >= Lớn hơn hoặc bằng Xuất dữ liệu: hàm fprintf, hàm disp  Hàm fprintf hay hàm disp được dùng để hiển thị lời giải thích hay lời ghi chú cho kết quả xuất ra trong Command Window Vd: Trong M-file, bạn gõ: fprintf(‘So SV lop 05VLUD-2 la:’) n Trong Command Window, bạn sẽ thấy xuất hiện: Sai cú pháp  Báo lỗi 3. Nhập dữ liệu và xuất dữ liệu a> Nhập dữ liệu: hàm input Vd: Trong M-file, bạn gõ: n=input(‘Nhap vao so SV lop 05VLUD-2:’) Trong Command Window, bạn sẽ được yêu cầu nhập vào giá trị của n (n là một số thực bất kỳ) Hoàn toàn tương tự với hàm disp.  Nếu trong bài lập trình của mình (tức là code bạn viết trong cửa sổ M-file) có những dòng bạn ghi chú hoặc giải thích thêm mà không VIC FAN
SWEET NOVEMBER 4 VERSION 2009 muốn Matlab thực hiện rồi báo lỗi, bạn hãy đặt fprintf(‘Loai yeu’) dấu % phía trước dòng đó. (Tác giả khuyên bạn elseif (diem==5)|(diem==6) nên sử dụng nhiều cái này vì tạo được kết cấu fprintf(‘Loai trung binh’) rõ ràng cho bài lập trình, đồng thời người đọc elseif(diem==7)|(diem==8) cũng thấy được “thiện ý”) fprintf(‘Loai kha’) Vd: Dấu * cho biết file có elseif(diem==9)|(diem==10) chỉnh sửa nhưng chưa lưu fprintf(‘Loai gioi’) else fprintf(‘So lieu khong hop le’) end (Ghi chú: để gõ dấu “|”, bạn nhấn đồng thời phím Shift và phím nằm ngay phía trên phím Enter) Nếu biểu thức logic của bạn chỉ có 2 giá trị thì trong cấu trúc lệnh if chỉ có if và else và Dòng bắt đầu bằng dấu %, có màu xanh lá cây end, không có else. Như ví dụ minh họa bên và Matlab không đọc những dòng này dưới.  Có sự khác nhau giữa dòng kết thúc bằng Vd: Nhập một số tự nhiên bất kỳ, lập trình để dấu “;” và dòng kết thúc không có gì cả. Nếu Matlab kiểm tra số đó có phải là số chẵn hay có dấu “;” thì kết quả của dòng đó không được không, xuất ra thông báo. xuất ra màn hình trong Command Window Để bắt đầu làm vd này, bạn cần làm quen với Vd: hàm lấy phần dư “mod”. Cấu trúc hàm mod Trong M-file, bạn gõ: như sau: a=1+2 mod(x,y) = n với x là số chia, y là số bị chia, n b=1+2; là phần dư khi lấy x chia y. Trong Command Window, bạn chỉ thấy: a= 3 (không hiển thị kết quả của b) 4. Lệnh điều kiện và vòng lặp a> Lệnh điều kiện: lệnh điều kiện hay dùng nhất trong Matlab là lệnh “if”. Lệnh “if: có cấu trúc như sau: if biểu thức logic 1 elseif biểu thức logic 2 Với vd đã cho, trong M-file bạn gõ đoạn code else sau: end so=input(‘Nhap vao so tu nhien can kiem Vd: Trong M-file, bạn gõ đoạn code sau, đây là tra:’) một chương trình xếp loại học lập dựa vào if (mod(so,2)==0) điểm số trung bình bạn nhập ban đầu. fprintf(‘So do la so chan’) diem=input(‘Hay nhap vao diem so trung else % Ở đây ngầm hiểu: mod(so,2)~=0 binh:’); fprintf(‘So do khong phai so chan’) if (diem>=1)&(diem<=4) end VIC FAN
SWEET NOVEMBER 5 VERSION 2009 đó trùng với ngày, tháng chuẩn đã gán thì hãy b> Vòng lặp: có 2 loại vòng lặp thường xuất ra câu “Chuc mung sinh nhat!”. (Thật là dùng trong Matlab: vòng lặp for và vòng lặp một ví dụ thú vị! Hè hè hè!) while. Chương trình của bạn như sau:  Vòng lặp for: được dùng khi sự lặp lại xảy namsinh=input(‘Nhap vao nam sinh cua ban:’); ra trong khoảng giới hạn xác định. Cấu trúc while (namsinh 2008) namsinh=input('Khong hop le. Vui long nhap lai vòng lặp for như sau: nam sinh cua ban:'); for a=a1:deltaa:a2 end % Neu ai do nhap nam sinh trong khoang 1930 den end % 1980 thi hop le, Matlab bo qua vong lap Giải thích: Cho biến a chạy từ giá trị ban đầu % while, nguoc lai neu % nam ngoai khoang 1930 % den 2008 thi vong lap while se bat ho phai nhap a1 đến giá trị cuối cùng a2 với bước nhảy là % lai. Tuong tu cho thang sinh va ngay sinh deltaa (a1, a2 và deltaa đều là các giá trị do bạn thangsinh=input(‘Nhap vao thang sinh cua ban:’); đưa vào). Nếu bạn không ghi deltaa thì Matlab while (thangsinh 12) sẽ ngầm hiểu deltaa=1 thangsinh=input('Khong hop le. Vui long nhap lai Vd: Giả sử bạn muốn xếp loại học tập của các thang sinh cua ban:'); end sinh viên lớp 05VLUD-2, bạn có thể dùng vòng ngaysinh=input(‘Nhap vao ngay sinh cua ban:’); lặp for như chương trình minh họa bên dưới: while (ngaysinh 31) ngaysinh=input('Khong hop le. Vui long nhap lai ngay sinh cua ban:'); end thangchuan=0; while (thangchuan 12) thangchuan=input('Hom nay la thang:'); end ngaychuan=0; while (ngaychuan 31) ngaychuan=input('Hom nay la ngay:'); end if (thangsinh==thangchuan)&(ngaysinh==ngaychuan) fprintf('CHUC MUNG SINH NHAT') end (Tất nhiên đây chỉ là ví dụ và ắt hẳn vẫn còn thiếu sót, như việc có tháng 31 ngày, có tháng 30 ngày, tháng 2 năm thường 29, năm nhuận  Vòng lặp while: được dùng khi sự lặp lại 28. Nếu bạn đọc thích có thể viết code thêm để xảy ra mà không xác định được khoảng giới xử lý vụ này) hạn. Nếu còn thỏa điều kiện lặp thì vòng lặp còn hoạt động, nếu điều kiện không còn thỏa 5. Một vài hàm thông dụng trong kỹ thuật thì thoát khỏi vòng lặp. Cấu trúc vòng lặp while mô phỏng bằng Matlab như sau: while biểu thức logic a> Hàm toán học sin: hàm sin end cos: hàm cos Vd: Viết một chương trình cho phép nhập tan: hàm tan ngày, tháng, năm sinh của người khác. Gán một abs: lấy giá trị tuyệt đối hoặc độ lớn của số ngày, tháng chuẩn bất kỳ, nếu sinh nhật của ai phức rem: hàm lấy phần dư sau khi chia VIC FAN
SWEET NOVEMBER 6 VERSION 2009 exp: hàm lũy thừa e log: logarit cơ số e log10: logarit cơ số 10 c> Hàm vẽ đồ thị Có khá nhiều hàm có thể dùng để vẽ đồ thị trong Matlab và tác giả cũng không am tường hết, nhưng trong đó hàm plot được ưu ái hơn cả. Hầu hết các bài mô phỏng, nếu có vẽ đồ thị, bạn đọc có thể dùng theo cấu trúc mẫu sau: plot(a,b) xlabel(‘Ten dai luong x (don vi)’) ylabel(‘Ten dai luong y (don vi)’) title(‘Ten do thi’) grid on Giải thích: plot(a,b): vẽ đồ thị biểu diễn sự thay đổi của đại lượng b theo đại lượng a (bắt buộc a và b phải là 2 mảng 1 chiều có cùng số phần tử!) xlabel: đặt tên cho trục hoành x của hệ tọa độ ylabel: đặt tên cho trục tung y của hệ tọa độ title(‘ ’): tên đồ thị grid on: chia lưới tọa độ (việc này giúp đồ thị của bạn dễ nhìn hơn) b> Hàm làm tròn Dưới đây là một ví dụ. round: làm tròn đến số nguyên gần nhất Vd: Vẽ đồ thị biểu diễn sự thay đổi điểm trung fix: làm tròn hướng về 0 bình học tập của một sinh viên qua 4 năm đại floor: làm tròn hướng xuống học theo bảng số liệu sau: ceil: làm tròn hướng lên Năm học Điểm trung bình học tập 2005 7.1 2006 5.9 VIC FAN
SWEET NOVEMBER 7 VERSION 2009 2007 6.7 2008 7.0 Code: namhoc=[2005 2006 2007 2008]; diemTB=[7.1 5.9 6.7 7.0]; plot(namhoc,diemTB) xlabel(‘Nam hoc’) ylabel(‘Diem trung binh hoc tap’) title(‘Do thi bieu dien su thay doi diem trung binh cua mot sinh vien qua 4 nam dai hoc’) grid on Khi nhấn menu Debug > Save and Run, một khung Figure sẽ nhảy ra cho bạn “nhân trạng” sau: Tại khung Save as type, chọn định dạng ảnh xuất ra là Bitmap files (*.bmp), đặt tên file rồi nhấn Save. (Bạn có thể thắc mắc sao không lưu với định dạng quen thuộc *.jpg ? Vì kinh nghiệm cho thấy ảnh xuất ra với đuôi này “mờ ảo” lắm, không đẹp đâu) Mặc định đường đồ thị của bạn sẽ có màu xanh dương (tiếc là trang giấy photo chỉ thể hiện được 2 màu trắng đen), bạn có thể chọn màu Nếu bạn có nhu cầu muốn xuất figure này ra khác tùy thích nếu theo cấu trúc: file ảnh (để chèn vào word báo cáo nộp thầy plot(a,b,‘r’) % duong do thi co mau do chẳng hạn),bạn vào menu File > Export plot(a,b,‘g’) % duong do thi mau xanh la cay Tương tự: b: màu xanh dương; c: màu lục lam; m: màu đỏ tươi; y: màu vàng; k: màu đen Ngoài ra, hàm plot còn hỗ trợ vẽ đồng thời nhiều đồ thị trên cùng một hệ tọa độ. Khi đó, cấu trúc của nó sẽ là: plot(a1,b1,a2,b2, ,an,bn) xlabel(‘Ten dai luong x (don vi)’) ylabel(‘Ten dai luong y (don vi)’) title(‘Ten do thi’) grid on Vd: Vẽ đồ thị biểu diễn sự thay đổi điểm trung bình học tập của 3 sinh viên qua 4 năm đại học theo bảng số liệu sau: Năm học Điểm trung bình học tập SV1 SV2 SV3 VIC FAN
SWEET NOVEMBER 8 VERSION 2009 2005 7.1 6.1 7.6 2006 5.9 5.7 7.4 2007 6.7 6.3 7.9 2008 7.0 5.9 8.3 Code: namhoc=[2005 2006 2007 2008]; diemTBSV1=[7.1 5.9 6.7 7.0]; diemTBSV2=[6.1 5.7 6.3 5.9]; diemTBSV3=[7.6 7.4 7.9 8.3]; plot(namhoc,diemTBSV1,namhoc,diemTB SV2,namhoc,diemTBSV3) xlabel(‘Nam hoc’)  rand(m,n): tạo ma trận m dòng n cột với các ylabel(‘Diem trung binh hoc tap’) phần tử có giá trị ngẫu nhiên nằm trong khoảng title(‘Do thi bieu dien su thay doi diem trung từ 0 đến 1 binh cua mot sinh vien qua 4 nam dai hoc’) Vd: grid on Figure: Mặc định, 3 đường sẽ có 3 màu khác nhau (lần  randint(m,n,[p,q]): tạo ma trận m dòng n cột lượt là xanh dương, đỏ, xanh lá), nếu bạn chưa với các phần tử có giá trị ngẫu nhiên là số ưng ý lắm có thể hiệu chỉnh màu theo cấu trúc: nguyên nằm trong khoảng [p,q] plot(namhoc,diemTBSV1,‘k’,namhoc,diem Vd: TBSV2,‘m’,namhoc,diemTBSV3,‘y’) d> Một số hàm khác  length(Y): xác định chiều dài của một mảng Y (tức là số phần tử có trong mảng Y) Vd: VIC FAN
SWEET NOVEMBER 9 VERSION 2009 6. Một vài thuật toán thông dụng trong kỹ thuật mô phỏng sử dụng Matlab a> Thuật toán tạo mảng 1 chiều Vd: Tạo mảng độ dày của màng Cách 1: Tạo trực tiếp doday=[345 786 890 299 102 999] (các phần tử cách nhau một khoảng trắng) Trong Command Window: doday= 345 786 890 299 102 999 Cách 2: Tạo bằng vòng lặp for N=input(‘Nhap vao so phan tu cua mang do day:’); for a=1:N doday(a)=input(‘Nhap do day:’); end doday Trong Command Window:  zeros(m,n): tạo ma trận m dòng n cột với các phần tử đều bằng 0.  eye(n): tạo ma trận đơn vi n dòng n cột  ones(m,n): tạo ma trận m dòng n cột với các phần tử đều bằng 1. b> Thuật toán tạo mảng 2 chiều (ma trận) Vd: Tạo ma trận m dòng n cột Cách 1: Tạo trực tiếp matran=[1 2 3;4 5 6] Trong Command Window: matran= 1 2 3 4 5 6 Cách 2: Tạo bằng 2 vòng lặp for lồng m=input(‘Nhap so dong cua ma tran:’); n=input(‘Nhap so cot cua ma tran:’); for a1=1:m VIC FAN
SWEET NOVEMBER 10 VERSION 2009 for a2=1:n d> Thuật toán sắp xếp mảng theo thứ matran(a1,a2)=input(‘Nhap phan tu:’) tự tăng dần (hoặc giảm dần) end Ý tưởng: Duyệt qua các phần tử có trong mảng, end so sánh giá trị của 2 phần tử cạnh nhau và hoán matran đổi vị trí của chúng. Trong Command Window: Vd: Tạo mảng gồm 10 số nguyên ngẫu nhiên nằm trong khoảng [1,50], sau đó sắp xếp mảng theo thứ tự tăng dần. % Tao mang ngau nhien mang=randint(1,10,[1,50]) % Sap xep mang theo thu tu tang dan for a1=1:length(mang)-1 for a2=a1+1:length(mang) if (mang(a1)>mang(a2)) % * tam=mang(a1); mang(a1)=mang(a2); mang(a2)=tam; end end end mang Trong Command Window: c> Thuật toán tổng cộng dồn Vd: Tính độ dày trung bình của màng doday=[157 890 456 228 456 761] % don vi nm N=length(doday); tong=0; % Tinh tong cua N phan tu do day (tong cong don) for a=1:N tong=tong+doday(a); end Nếu muốn sắp xếp mảng theo thứ tự giảm dần, tong % Tinh do day trung binh bạn chỉ việc thay đổi dòng % * ở trên thành: dodaytrungbinh=tong/N if (mang(a1)<mang(a2)) Trong Command Window: Khi đó: VIC FAN
SWEET NOVEMBER 11 VERSION 2009 e> Thuật toán chuyển số thập phân for a=1:sobit sang số nhị phân và ngược lại sonhiphan(a)=input('Nhap vao tung bit % Thuat toan chuyen tu so thap phan sang cua so nhi phan:'); so nhi phan end sothapphan=input('Nhap vao so thap phan sonhiphan can chuyen doi:') tong=0; chay=1; dem=sobit-1; while round(sothapphan/2)~=0 for b=1:sobit sodu(chay)=rem(sothapphan,2); tong=tong+sonhiphan(b)*2^dem; chay=chay+1; dem=dem-1; sothapphan=floor(sothapphan/2); end end sothapphan=tong % Nghich dao mang sodu ta duoc mang so Trong Command Window: nhi phan tam=0; for a=1:length(sodu) sonhiphan(length(sodu)-tam)=sodu(a); tam=tam+1; end sonhiphan Trong Command Window: % Thuat toan chuyen tu so nhi phan sang so thap phan sobit=input('Nhap vao so bit cua so nhi phan can chuyen:'); VIC FAN
SWEET NOVEMBER 12 VERSION 2009 7. Một vài lưu ý trong sử dụng Matlab 8. Một vài ví dụ giúp thực hành các thuật toán  Matlab phân biệt ký tự viết hoa và ký tự  Viết chương trình giúp giải phương trình viết thường. bậc 2  Phần đầu mỗi bài lập trình, bạn nên có 2 clc dòng lệnh sau: clear all clc a=input('Nhap vao he so a cua phuong clear all trinh bac 2:'); clc là lệnh xóa màn hình trong Command b=input('Nhap vao he so b cua phuong Window, tạo giao diện “thoáng mát” và trinh bac 2:'); “tươi mới” cho mỗi lần chạy. c=input('Nhap vao he so c cua phuong clear all là lệnh xóa tất cả các biến đã gán trinh bac 2:'); trong chương trình, đảm bảo kết quả của lần delta=b^2-4*a*c; chạy trước không gây ảnh hưởng đến lần if (delta 0 (clear, a, b đều cách nhau một khoảng fprintf('Nghiem phuong trinh la:') trắng) x1=(-b+sqrt(delta))/2  Hiển thị kết quả trong Matlab: nếu kết x2=(-b-sqrt(delta))/2 quả của bạn là 0.001 thì Matlab vẫn để end nguyên là 0.001. Nhưng nếu bắt đầu lấn sau 10-4 thì Matlab bắt đầu “giở trò” a=0.0001 sẽ được viết thành a=1.0000e-004 Vậy là e-004 chính là nhân 10-4. Nhớ ý này nha bạn!  Đặt tên file: Khi bạn save lại đoạn chương trình hay bài lập trình của mình, tên file được đặt không được có dấu (hay ký hiệu đặc biệt), nếu không khi bạn Run file này, Matlab sẽ có những báo lỗi khó hiểu. Và thêm 1 điều nữa, để tránh “xung đột” đáng tiếc có thể xảy ra, trong quá trình lập trình bằng Matlab bạn hãy tắt Unikey hay Vietkey đi (chúng hiếm khi hòa hợp nhau lắm).  “Chuyển nhà” cho file .m: Khi bạn copy, cut, paste một file Matlab (file có đuôi .m) từ nơi này sang nơi khác, rất lưu ý là muốn Matlab chạy được file này trơn tru thì bạn phải đặt file này vào folder work của Matlab (folder work thường nằm trong folder MATLAB6p5, folder MATLAB6p5 có thể nằm trong Program Files của ổ đĩa mà bạn cài Matlab). Hết sức lưu ý! VIC FAN
SWEET NOVEMBER 13 VERSION 2009  Viết chương trình tạo mảng n phần tử có số clear all hạng thứ i là tổng của hai số hạng i – 1 và i – 2. n=input('Nhap vao so dong cua ma tran:'); clc m=input('Nhap vao so cot cua ma tran:'); clear all l=1; n=input('Nhap vao so phan tu co trong % Nhap ma tran A va chuyen ma tran A mang:'); sang mang B a=input('Nhap vao phan tu thu nhat cua for r=1:n mang:'); for s=1:m mang(1)=a; A(r,s)=input('Nhap vao cac phan tu b=input('Nhap vao phan tu thu hai cua trong ma tran:'); mang:'); B(l)=A(r,s); mang(2)=b; l=l+1; i=3; end % Cach 1: Dung vong lap for end for bien=3:n % "bien" chay tu 3 den n buoc nhay la 1 fprintf('Ma tran A truoc khi sap xep:') mang(bien)=mang(bien-2)+mang(bien-1); A end % Sap xep cac phan tu trong mang B theo mang thu tu tang dan % Cach 2: Dung vong lap while for w=1:m*n-1 bien=3; for v=w+1:m*n while bien B(v)) % * mang(bien)=mang(bien-2)+mang(bien-1); tam=B(w); bien=bien+1; B(w)=B(v); end B(v)=tam; mang end end end % Chuyen mang B sau khi sap xep thanh ma tran A h=1; for g=1:n for f=1:m A(g,f)=B(h); h=h+1; end end fprintf('Ma tran A sau khi sap xep:')  Viết chương trình sắp xếp ma trận sao cho A các phần tử thay đổi theo thứ tự tăng dần. Ví (Nếu sắp xếp theo thứ tự giảm dần thì ở dòng dụ: % *, câu lệnh sẽ là: 89 27 0 0 1 5 if (B(w)<B(v)) ) A = 1 102 5 → A = 6 11 27 6 11 34 34 89 102 Ý tưởng: Chuyển ma trận sang mảng; Sắp xếp mảng theo thứ tự tăng dần; Chuyển mảng thành ma trận. clc VIC FAN
SWEET NOVEMBER 14 VERSION 2009 fprintf('Ma tran sau khi chuyen doi:'); b  Viết chương trình chuyển đổi cột thành dòng trong ma trận. Ví dụ: 2 16 23 2 6 7 2 6 1 15 a = 16 1 8 22 → a = 7 8 6 23 15 6 5 2 22 5 clc clear all n=input('Nhap vao so dong cua ma tran:'); m=input('Nhap vao so cot cua ma tran:'); for r=1:n 9. Một số bài tập mẫu for s=1:m Bài tập 1: Tạo ma trận. Bạn hãy lần lượt thực a(r,s)=input('Nhap vao cac phan tu hiện các yêu cầu sau: trong ma tran:'); a) Tạo một ma trận A là ma trận m dòng, n cột b(s,r)=a(r,s); với các phần tử là số nhập vào bất kỳ. end clc end clear all fprintf('Ma tran truoc khi chuyen doi:') m=input('Nhap vao so dong cua ma tran:'); a n=input('Nhap vao so cot cua ma tran:'); % Nhap vao cac phan tu cua ma tran VIC FAN
SWEET NOVEMBER 15 VERSION 2009 for a=1:m for b=1:n A(a,b)=input('Nhap phan tu:'); end end A Ví dụ: d) Tạo một ma trận D là ma trận m dòng n cột với các phần tử là số thập phân ngẫu nhiên bất kỳ nằm từ 0 đến 1. clc b) Tạo một ma trận B là ma trận vuông n với clear all các phần tử là số nhập vào bất kỳ. m=input('Nhap vao so dong cua ma tran:'); clc n=input('Nhap vao so cot cua ma tran:'); clear all D=rand(m,n) n=input('Nhap vao so cot hay dong cua ma Ví dụ: tran:'); % Nhap vao cac phan tu cua ma tran for c=1:n for d=1:n B(c,d)=input('Nhap phan tu:'); end end B Ví dụ: e) Tạo một ma trận E là ma trận m dòng n cột với các phần tử là số thập phân ngẫu nhiên bất kỳ nằm từ 10 đến 100. clc clear all c) Tạo một ma trận C là ma trận m dòng n cột m=input('Nhap vao so dong cua ma tran:'); với các phần tử là số nguyên ngẫu nhiên bất kỳ n=input('Nhap vao so cot cua ma tran:'); nằm từ 10 đến 100. E=rand(m,n)+randint(m,n,[10,100-1]) clc Ví dụ: clear all m=input('Nhap vao so dong cua ma tran:'); n=input('Nhap vao so cot cua ma tran:'); C=randint(m,n,[10,100]) Ví dụ: f) Tạo ma trận F là ma trận m dòng n cột có các phần tử ngẫu nhiên đều là số chẵn và ma trận G VIC FAN
SWEET NOVEMBER 16 VERSION 2009 là ma trận m dòng n cột có các phần tử ngẫu end nhiên đều là số lẻ. (Các phần tử nằm trong end khoảng từ 10 đến 100) H clc Ví dụ: clear all m=input('Nhap vao so dong cua ma tran:'); n=input('Nhap vao so cot cua ma tran:'); F=2*randint(m,n,[10/2,100/2]) G=2*randint(m,n,[10/2,100/2-1])+1 Ví dụ: Bài tập 2: Thao tác trên ma trận. Bạn hãy lần lượt thực hiện các yêu cầu sau: a) Tạo ma trận A là một mảng 1 chiều m phần tử. Nhập vào các phần tử bất kỳ. Tìm phần tử lớn nhất, phần tử nhỏ nhất, tính giá trị trung bình của các phần tử. clc clear all m=input('Nhap vao so phan tu cua mang A:'); % Cach 1: Dung ham co san (De! Don gian!) for a=1:m g) Tạo một ma trận H là ma trận m dòng n cột A(a)=input('Nhap phan tu:'); với dòng chứa phần tử chẵn và dòng chứa phần end tử lẻ xen kẽ với nhau ((Các phần tử nằm trong fprintf('So lon nhat trong mang A:') khoảng từ 10 đến 100). Ví dụ: sln=max(A) fprintf('So nho nhat trong mang A:'); H = snn=min(A) Dòng chứa phần tử chẵn fprintf('Gia tri trung binh cua cac phan tu:'); trungbinh=mean(A) Dòng chứa phần tử lẻ % Cach 2: Dung lap trinh (Kho nhung tri tue hon!) sln=0; clc snn=10000; % Gan snn ban dau la mot so rat clear all lon m=input('Nhap vao so dong cua ma tran:'); tong=0; n=input('Nhap vao so cot cua ma tran:'); for a=1:m % Nhap vao cac phan tu cua ma tran A(a)=input('Nhap phan tu:'); for i=1:m if (sln A(a)) % Tim so nho nhat H(i,j)=2*randint(1,1,[10/2,100/2-1])+1; snn=A(a); else % i la so le, tuc dong chua phan tu end chẵn tong=tong+A(a); H(i,j)=2*randint(1,1,[10/2,100/2]); end end VIC FAN
SWEET NOVEMBER 17 VERSION 2009 sln end snn B trungbinh=tong/m fprintf('So lon nhat trong mang B la:') sln=max(max(B)) fprintf('So nho nhat trong mang B la:') snn=min(min(B)) % Cach 2: Dung lap trinh sln=0; snn=10000; for a=1:n for b=1:n B(a,b)=input('Nhap phan tu:'); if (sln B(a,b)) % Tim so nho nhat không ghi nhận được sự tồn tại của mảng A (A snn=B(a,b); bây giờ chỉ là một số): end sln=0; end snn=10000; % Gan snn ban dau la mot so rat end lon sln tong=0; snn for a=1:m Ví dụ: A=input('Nhap phan tu:'); if (sln A) % Tim so nho nhat snn=A; end tong=tong+A; end sln snn trungbinh=tong/m b) Tạo ma trận B là ma trận vuông n. b1> Tìm phần tử lớn nhất, phần tử nhỏ nhất trong ma trận. clc b2> Tạo f là một số ngẫu nhiên từ 1 đến n, clear all hãy xuất ra dòng f, cột f của ma trận. n=input('Nhap vao so dong hay so cot cua % Tao so ngau nhien f ma tran:'); f=randint(1,1,[1,n]) % Cach 1: Dung ham % Xuat ra dong f cua ma tran B for a=1:n fprintf('Dong f cua ma tran B la:') for b=1:n dongf=B(f,:) % B(f,:) nghia la tat ca cac cot B(a,b)=input('Nhap phan tu:'); nam o dong f end % Xuat ra cot f cua ma tran B VIC FAN
SWEET NOVEMBER 18 VERSION 2009 fprintf('Cot f cua ma tran B la:') cotf=B(:,f) % B(:,f) nghia la tat ca cac dong nam o cot f Ví dụ: b3> Tạo s là một số ngẫu nhiên khác từ 1 đến n (s ≠ f). Hãy đổi chỗ qua lại giữa dòng f và dòng s, giữa cột f và cột s. s=randint(1,1,[1,n]); while (s==f) % Vong lap bat buoc s phai khac f s=randint(1,1,[1,n]); end s dongs=B(s,:); cots=B(:,s); % Doi dong f va dong s fprintf('Doi dong f va dong s:') b4> Chèn thêm dòng thứ n + 1 vào ma trận B(f,:)=dongs; B ở b3 có các phần tử là tổng của các phần tử B(s,:)=dongf; tương ứng giữa dòng thứ nhất và dòng thứ n. B dong1=B(1,:) % Doi cot f va cot s dongn=B(n,:) fprintf('Doi tiep cot f va cot s:') tong=dong1+dongn cotf=B(:,f); fprintf('Ma tran B sau khi them dong n + 1:') cots=B(:,s); B(n+1,:)=tong B(:,f)=cots; Ví dụ: B(:,s)=cotf; B Ví dụ: VIC FAN
SWEET NOVEMBER 19 VERSION 2009 c) Tạo ma trận C là ma trận vuông n. c1> Hãy xuất ra các phần tử trên đường chéo chính, đường chéo phụ của ma trận. clc clear all b5> Cho biết có bao nhiêu số chẵn, số lẻ n=input('Nhap vao so dong hay so cot cua ma trong ma trận mới ở b4. tran C:'); dem1=0; % TIM DUONG CHEO CHINH dem2=0; % Cach 1: Dung ham for c=1:n for a=1:n for d=1:n for b=1:n if (mod(B(c,d),2)==0) % Tim so chan C(a,b)=input('Nhap phan tu:'); dem1=dem1+1; end end end duongcheochinh=diag(C) end % Cach 2: Dung lap trinh end dem=1; fprintf('So phan tu la so chan:') for a=1:n dem1 for b=1:n fprintf('So phan tu la so le:') C(a,b)=input('Nhap phan tu:'); n*n-dem1 end Ví dụ: duongcheochinh(dem)=C(a,a); dem=dem+1; end C duongcheochinh Ví dụ: VIC FAN
SWEET NOVEMBER 20 VERSION 2009 c2> Tìm phần tử lớn nhất, nhỏ nhất trên đường chéo chính, đường chéo phụ của ma trận. Câu này hoàn toàn tương tự như câu a, bạn có thể dùng hàm min, max rất nhanh và đơn giản. Ở đây giới thiệu đoạn code trong trường hợp phải lập trình: % DUONG CHEO CHINH sln1=0; snn1=10000; for a1=1:length(duongcheochinh) if (sln1 duongcheochinh(a1)) dem2=0; snn1=duongcheochinh(a1); for a=1:n end for b=1:n end sln1 C(a,b)=input('Nhap phan tu:'); snn1 end % DUONG CHEO PHU duongcheophu(dem1)=C(a,n-dem2); sln2=0; dem1=dem1+1; snn2=10000; dem2=dem2+1; for a2=1:length(duongcheophu) end if (sln2 duongcheophu(a2)) snn2=duongcheophu(a2); % TIM DUONG CHEO PHU end dem1=1; end dem2=0; sln2 for a=1:n snn2 duongcheophu(dem1)=C(a,n-dem2); Ví dụ: dem1=dem1+1; dem2=dem2+1; end C duongcheophu Ví dụ: c3> Sắp xếp các phần tử trên đường chéo chính theo thứ tự tăng dần. Xuất ra ma trận C sau khi sắp xếp. % DUONG CHEO CHINH % Sap xep cac phan tu tang dan for g=1:length(duongcheochinh)-1 for h=g+1:length(duongcheochinh) VIC FAN
SWEET NOVEMBER 21 VERSION 2009 if (duongcheochinh(g)>duongcheochinh(h)) tam=duongcheochinh(g); duongcheochinh(g)=duongcheochinh(h); duongcheochinh(h)=tam; end end end duongcheochinh % Ma tran C sau khi sap xep fprintf('Ma tran C sau khi sap xep:'); c5> Giả sử n là số chẵn, hãy xuất ra ma trận for k=1:length(duongcheochinh) C(k,k)=duongcheochinh(k); C có các phần tử trên đường chéo chính đều end bằng 0 và các phần tử trên đường chéo phụ đều C bằng 1. Ví dụ: % Dua phan tu 0 vao duong cheo chinh for w=1:n C(w,w)=0; end % Dua phan tu 1 vao duong cheo phu t=0; for v=1:n C(v,n-t)=1; t=t+1; end c4> Từ ma trận C ban đầu (khi chưa làm C c3), sắp xếp các phần tử trên đường chéo phụ theo thứ tự giảm dần. Xuất ra ma trận C sau khi sắp xếp. % DUONG CHEO PHU % Sap xep cac phan tu giam dan for g=1:length(duongcheophu)-1 for h=g+1:length(duongcheophu) if (duongcheophu(g)<duongcheophu(h)) d) Nhập vào một ma trận vuông, xuất ra phần tam=duongcheophu(g); tử lớn nhất trên từng dòng, sau đó đổi chỗ các duongcheophu(g)=duongcheophu(h); phần tử này vào đường chéo chính của ma trận. duongcheophu(h)=tam; clc end clear all end n=input('Nhap so dong hay so cot cua ma end tran:'); duongcheophu m=n; % Ma tran C sau khi sap xep for i=1:n fprintf('Ma tran C sau khi sap xep:'); for j=1:m q=0; A(i,j)=input('Nhap phan tu:'); for p=1:length(duongcheophu) end C(p,n-q)=duongcheophu(p); end q=q+1; A end % Tim phan tu lon nhat tren tung dong C dem=1; Ví dụ: for i=1:n VIC FAN
SWEET NOVEMBER 22 VERSION 2009 maxdong(i)=max(A(i,:)); % Tao ma tran A for j=1:m for a=1:m if (A(i,j)==maxdong(i)) for b=1:n vitri(dem)=j; A(a,b)=input('Nhap phan tu ma tran A:'); dem=dem+1; end end end end A end % Tao ma tran B maxdong for a=1:m vitri for b=1:n % Dua phan tu lon nhat tung dong ve duong B(a,b)=input('Nhap phan tu ma tran B:'); cheo chinh end B=A; end i=1; B for j=1:n Ví dụ: A(i,i)=maxdong(i); A(i,vitri(j))=B(i,j); i=i+1; end A Ví dụ: a1> Hãy tính ma trận C là tổng của ma trận A và ma trận B. fprintf('Ma tran C:') C=A+B Ví dụ: a2> Hãy tính ma trận D là hiệu của ma trận A và ma trận B. fprintf('Ma tran D:') Bài tập 3: Các phép tính trên ma trận D1=A-B a) Tạo ma trận A và ma trận B đều là ma trận m D2=B-A dòng, n cột. Ví dụ: clc clear all m=input('Nhap vao so dong cua ma tran A:'); n=input('Nhap vao so cot cua ma tran A:'); VIC FAN
SWEET NOVEMBER 23 VERSION 2009 a3> Hãy tính ma trận F là tích của các phần tử tương ứng của ma trận A và ma trận B. fprintf('Ma tran tich phan tu F:') for a=1:m for b=1:n F(a,b)=A(a,b)*B(a,b); end c) Thêm dòng hay cột vào các ma trận A, B và end G để chúng trở thành ma trận vuông. Tính ma F trận I là tích của 3 ma trận A, B, G (A*B*G; B*G*A; G*A*B; A*G*B; G*B*A; B*A*G;) Bài tập 4: (Bài tập bổ sung) a) Nhập vào một số tự nhiên bất kỳ. Viết b) Tạo ma trận G là ma trận n dòng, m cột. Hãy chương trình kiểm tra số đó có phải là số tính ma trận H là tích của ma trận A và ma trận nguyên tố hay không. G (A*G và G*A). clc clear all % Tao ma tran G N=input('Nhap vao so tu nhien bat ky:') for c=1:n kiemtra=1; for d=1:m for a=2:N-1 G(c,d)=input('Nhap phan tu ma tran if (mod(N,a)==0) G:'); kiemtra=0; end end end end G if (kiemtra==0) % Tinh tich ma tran A va G fprintf('So nhap vao khong phai so nguyen H1=A*G to') H2=G*A else fprintf('So nhap vao la so nguyen to') Ví dụ: end b) Nhập ma trận A gồm các số tự nhiên. Hãy xuất ra các số nguyên tố có trong ma trận A. Cho biết có bao nhiêu số nguyên tố có trong ma trận A. (Ghi chú: Số nguyên tố là số chỉ chia hết cho 1 và chính nó) clc VIC FAN
SWEET NOVEMBER 24 VERSION 2009 clear all m=input('Nhap vao so dong cua ma tran:'); n=input('Nhap vao so cot cua ma tran:'); for a=1:m for b=1:n A(a,b)=input('Nhap phan tu:'); end end B=A; % Kiem tra phan tu nao la so nguyen to for a=1:m for b=1:n for s=2:(A(a,b)-1) if (mod(A(a,b),s)==0) B(a,b)=0; end end end end % Xuat cac phan tu la so nguyen to fprintf('Cac so nguyen to trong ma tran A la:') B. HƯỚNG DẪN THỰC HÀNH CÁC BÀI dem=1; MÔ PHỎNG for a=1:m BÀI LẬP TRÌNH 1 for b=1:n XÁC ĐỊNH CHIẾT SUẤT VÀ ĐỘ DÀY if (B(a,b)~=0) CỦA MÀNG TỪ PHỔ TRUYỀN QUA snt(dem)=B(a,b); Cho bảng số liệu độ truyền qua T tương ứng dem=dem+1; với bước sóng λ(nm) và phổ truyền qua của end một mẩu màng đa lớp như sau: end end snt % Xoa nhung so nguyen to trung nhau for c=1:length(snt)-1 for d=c+1:length(snt) if (snt(c)==snt(d)) snt(d)=0; end end end dem=1; for e=1:length(snt) if (snt(e)~=0) snt_new(dem)=snt(e); dem=dem+1; end end snt_new % Dem so nguyen to fprintf('So so nguyen to co ma tran A:') length(snt_new) Ví dụ: VIC FAN
SWEET NOVEMBER 25 VERSION 2009 Các bạn thực hiện lần lượt các yêu cầu sau: những vân sáng tối xen kẽ nhau, những vân 1) Nội suy các giá trị độ truyền qua T Mi và Tmi sáng tương ứng với những giá trị T cực đại và ứng với tất cả các bước sóng λ i tại các giá trị những vân tối tương ứng với những giá trị T cực đại và cực tiểu của phổ. cực tiểu. Từ phổ truyền qua T(λ), người ta có 2) Xác định chiết suất n 1 của màng theo công thể xác định được độ dày và chiết suất của thức (14) ứng với tất cả các cực trị (Cho chiết màng theo phương pháp Swanepoel. Trong suất đế S = 1.52) phương pháp này, sẽ có những giá trị T ảo 1/ 2 không có trên phổ được nội suy để hỗ trợ cho CT (14): n N N 2 S 2 việc tính toán. T T S 2 1 Bước đầu tiên bạn nhập trong M-file mấy Trong đó: N 2S M m dòng sau, chính là bảng số liệu thầy cho các TMTm 2 bạn: 3) Xác định độ dày di của màng theo công thức Tmax0 (26) ứng với tất cả các cực trị.   CT (26): d 1 2 2 1n2 2n1 lamdamax Trong đó: n1 và n2 là chiết suất ở hai giá trị bước sóng liên tiếp. 4) Tính độ dày trung bình của của các di vừa tìm, ta tạm gọi là độ dày d1. 5) Dùng n1 và d1 để xác định bậc m của các cực trị theo công thức (5) lamdamin Tmin0 CT (5): 2nd = mλ % So lieu tu pho truyen qua 6) Dùng lại n 1 và m để tính lại độ dày d 2 dựa % Mang Tmax0 ung voi cac dinh cuc dai tren pho theo công thức (5) Tmax0=[87.375 86.935 86.911 86.758 86.718 86.363 86.145 85.056 84.915 84.637 83.557 81.11 7) Từ câu 6, tính lại chiết suất n2. 79.309 74.851 69.294]; Sau khi viết xong chương trình, chạy ổn định % Mang lamdamax tuong ung voi cac cuc dai trong và in ra bảng kết quả: mang Tmax0 lamdamax=[692 604 645 749 569 537 510 813 484 892 989 462 443 426 410]; % Mang Tmin0 ung voi cac dinh cuc tieu tren pho Hướng dẫn: Tmin0=[3.2298*10^-3 47.336 56.142 62.778 66.675 Đây là bài toán thuận và tương đối đơn giản 71.055 71.725 72.47 72.902 73.8 74.349 74.359 nhất trong loạt bài mô phỏng quang học (nói 75.127 75.161 75.766]; % Mang lamdamin tuong ung voi cac cuc tieu trong vậy chứ đối với những ai mới bắt đầu đều phải mang Tmin0 “lè lưỡi”). Trước hết chúng ta sẽ tìm hiểu ý lamdamin=[304 386 401 416 433 452 496 472 939 nghĩa của phổ truyền qua, nó đơn giản thế này: 851 719 523 552 781 667]; I  Tuy nhiên, chúng ta không cần phải lấy hết T  d I  toàn bộ vùng phổ, chúng ta chỉ lấy phần phổ o truyền qua có dáng đẹp nhất (kết quả sẽ chính Trong đó I (λ) là cường độ ánh sáng tới o xác hơn), phần mà nó uốn lượn hình sin đều mẫu; I (λ) là cường độ ánh sáng sau khi đi qua d đặn, thường là vùng bước sóng khả kiến. Do mẫu có độ dày d. vậy, bạn thêm phần chương trình bên dưới để Phổ truyền qua T(λ) có hình dạng uốn lượn lọc lại những phần tử đúng. như trong bảng số liệu bên trái. Tại sao như % Chon vung mo phong la vung kha kien, tuc buoc vậy? À, thì ra là do hiện tượng giao thoa ánh song tu 400 den 750nm sáng khi ánh sáng đi qua bản mỏng (ở đây là t1=1; màng mỏng của chúng ta). Sự giao thoa tạo nên for k1=1:length(lamdamax) VIC FAN
SWEET NOVEMBER 26 VERSION 2009 if (lamdamax(k1)>400)&(lamdamax(k1) 400)&(lamdamin(k2) lamdamax(f2)) clear lamdamin Tmin0 F=Tlamda_max(:,f1); lamdamin=lamdamin_1 Tlamda_max(:,f1)=Tlamda_max(:,f2); Tmin0=Tmin0_1 Tlamda_max(:,f2)=F; A, hình như 2 mảng bước sóng của chúng lamdamax=Tlamda_max(1,:); ta không được trật tự lắm, chúng cần được sắp end xếp tăng dần (hoặc giảm dần) để có thể dễ dàng end end thao tác. Điều đáng nói ở đây là khi bạn sắp lamdamax xếp mảng lamdamax thì đồng thời cũng phải Tmax0=Tlamda_max(2,:) sắp xếp mảng Tmax0 tương ứng. Tương tự, khi % Mang lamdamin va Tmin0 bạn sắp xếp mảng lamdamin thì đồng thời cũng Tlamda_min=[lamdamin;Tmin0] phải sắp xếp mảng Tmin0 tương ứng. Việc lập for g1=1:length(lamdamin)-1 for g2=g1+1:length(lamdamin) trình vụ này “hơi khó”, bạn có thể lập trình thế if (lamdamin(g1)>lamdamin(g2)) này (có 3 cách cho các bạn lựa chọn): G=Tlamda_min(:,g1); Cách 1: Tlamda_min(:,g1)=Tlamda_min(:,g2); % Sap xep cac mang buoc song theo thu tu tang Tlamda_min(:,g2)=G; dan lamdamin=Tlamda_min(1,:); % Mang lamdamax va Tmax0 end for a1=1:length(lamdamax)-1 end for b1=a1+1:length(lamdamax) end if (lamdamax(a1)>lamdamax(b1)) lamdamin tam1=lamdamax(a1); Tmin0=Tlamda_min(2,:) lamdamax(a1)=lamdamax(b1); Cách 3: Dùng hàm chuyên sắp xếp mảng – hàm lamdamax(b1)=tam1; sort (Dễ nhất, đơn giản nhất và ngắn gọn nhất) tam2=Tmax0(a1); Tmax0(a1)=Tmax0(b1); % Sap xep cac mang buoc song theo thu tu tang Tmax0(b1)=tam2; dan end % Mang lamdamax va Tmax0 end lamda_max1=lamdamax end lamdamax=sort(lamdamax) lamdamax Tmax0 for thu1=1:length(lamdamax) % Mang lamdamin va Tmin0 for thu2=1:length(lamdamax_1) for a2=1:length(lamdamin)-1 if (lamdamax(thu1)==lamdamax_1(thu2)) for b2=a2+1:length(lamdamin) tam=Tmax0(thu1)=Tmax0(thu2); if (lamdamin(a2)>lamdamin(b2)) end tam3=lamdamin(a2); VIC FAN
SWEET NOVEMBER 27 VERSION 2009 end for v=1:length(lamdamax) end Tmin1(v)=interp1(lamdamin,Tmin0, Tmax0 lamdamax(v)); 4 mảng trên sau khi sắp xếp bây giờ thế này: end lamdamax = Tmin1 410 426 443 462 484 510 537 569 604 Bạn sẽ được mảng Tmin1 như sau: 645 692 749 Tmin1 = Tmax0 = 60.1236 65.0704 68.9803 71.7625 69.2940 74.8510 79.3090 81.1100 84.9150 72.0975 73.0908 74.7298 75.2215 86.1450 86.3630 86.7180 86.9350 86.9110 75.4159 75.6438 75.0848 NaN 87.3750 86.7580 Á, NaN là gì? Đừng nao núng, hãy tiếp tục lamdamin = làm tương tự với Tmax1 401 416 433 452 472 496 523 552 667 % Tim Tmax1 tai cac gia tri buoc song 719 lamdamin đa biet Tmin0. Noi suy duoc thuc Tmin0 = hien dua tren cac cap gia tri 56.1420 62.7780 66.6750 71.0550 72.4700 71.7250 74.3590 75.1270 75.7660 74.3490 (lamdamax,Tmax0) da biet tren pho Bước tiếp theo chúng ta cần làm, cũng là for v=1:length(lamdamin) Tmax1(v)=interp1(lamdamax,Tmax0, khâu quan trọng nhất, khó nhất trong bài lập lamdamin(v)); trình này, chính là nội suy các giá trị T ảo end không có trên phổ. Tại giá trị bước sóng nào đã Tmax1 có T cực đại (Tmax0) thì ta phải tìm ta T cực Bạn sẽ được mảng Tmax1 như sau: tiểu (Tmin1) tại bước sóng đó. Tương tự, tại Tmax1 = giá trị bước sóng nào đã có T cực tiểu (Tmin0) NaN 71.3779 76.6866 80.1621 82.8395 thì ta phải nội suy tìm T cực đại (Tmax1). 85.4827 86.2500 86.5294 87.1282 87.0827 Thế NaN là gì? Tại sao chúng xuất hiện một cách bất thường như vậy? NaN chẳng qua chỉ là giá trị không xác định thôi. Nội suy interp1 có giới hạn của nó. Nó không thể nội suy được D F tất cả, đặc biệt là những giá trị nằm ở biên đầu E B B và biên cuối. Công việc của chúng ta tiếp theo A là “khai trừ” những phần tử đáng ghét này. Để B C C B B tạo sự phong phú và đa dạng cho bài lập trình (tránh người trên nghi ngờ), ở đây tác giả cũng giới thiệu các bạn 2 cách. Theo quy tắc, giá trị tại điểm A sẽ được nội % Khu bo NaN suy từ 2 giá trị lân cận B và C (Ở đây giá trị % Cach 1 (đây là cách thủ công, đơn giản, chỉ điểm A thuộc mảng Tmin1). Tương tự, giá trị áp dụng được trong bài này, khi NaN nằm ở tại điểm D sẽ được nội suy từ 2 giá trị lân cận E đầu mảng Tmin1 và cuối mảng Tmax1) và F (Ở đây giá trị điểm D thuộc mảng Tmax1). Tmin1(length(Tmin1))=[]; % [] la ky hieu xoa Chúng ta sẽ dùng hàm interp1 để tiến hành phan tu nội suy. Bạn đọc có thể tham khảo đoạn Tmax1(1)=[]; chương trình sau: Tmin1 % Cau 1 Tmax1 % Tim Tmin1 tai cac gia tri buoc song Bạn sẽ được mảng Tmin1 và Tmax1 sau khi lamdamax da biet Tmax0. Noi suy duoc thuc khử NaN như sau: hien dua tren cac cap gia tri (lamdamin,Tmin0) Tmin1 = da biet tren pho VIC FAN
SWEET NOVEMBER 28 VERSION 2009 60.1236 65.0704 68.9803 71.7625 Tmax0=mang1b 72.0975 73.0908 74.7298 75.2215 % Khu bo NaN trong mang Tmax1 75.4159 75.6438 75.0848 B=isfinite(Tmax1) Tmax1 = dem2=1; 71.3779 76.6866 80.1621 82.8395 for bien2=1:length(Tmax1) 85.4827 86.2500 86.5294 87.1282 87.0827 if (B(bien2)==1) mang2(dem2)=Tmax1(bien2); Ồ, hay quá, mất rồi. Tuy nhiên, vấn đề mang2a(dem2)=lamdamin(bien2); không đơn giản là việc bạn bỏ phần tử này, giữ mang2b(dem2)=Tmin0(bien2); phần tử kia. Nếu bạn đã xóa đi 1 phần tử trong dem2=dem2+1; mảng Tmin1 thì bạn cũng nên xóa đi phần tử ở end vị trí tương ứng trong mảng lamdamax end lamdamax(length(lamdamax))=[]; clear Tmax1 lamdamin Tmin0 Tương tự với mảng lamdamin vốn tương Tmax1=mang2 ứng với Tmax1 lamdamin=mang2a lamdamin(1)=[]; Tmin0=mang2b Không những thế, bạn cũng phải xóa luôn Ngoài hàm isfinite, bạn đọc cũng có thể phần tử tương ứng trong mảng Tmin0. Tại sao? dùng hàm “chuyện trị” isnan với cách tác động Vì Tmin0 liên hệ với lamdamin rất mật thiết như sau: (chúng là cặp “bài trùng” đi cùng nhau trên phổ isnan(NaN)=1 mà). Việc xóa bỏ này đảm bảo tính tương thích. isnan(Số bất kỳ)=0 Tmin0(1)=[]; Như vậy bạn đã có thêm cách thứ 3 rồi đấy. Tương tự với mảng Tmax0 vốn tương thích % Cach 3 với lamdamax % Khu bo NaN trong mang Tmin1 Tmax0(length(lamdamax))=[]; A=isnan(Tmin1) % Cach 2 (cách này khoa học hơn, tổng quát dem1=1; hơn và “đáng nể hơn”, có thể áp dụng cho for bien1=1:length(Tmin1) nhiều bài khác, khai trừ NaN ở mọi vị trí bất kỳ if (A(bien1)==0) mang1(dem1)=Tmin1(bien1); trong mảng. Đó là việc dùng hàm isfinite, một mang1a(dem1)=lamdamax(bien1); anh chàng chuyên trị các nàng NaN (không xác mang1b(dem1)=Tmax0(bien1); định) và inf (vô cùng) đỏng đảnh, khó chịu. Tác dem1=dem1+1; động của hàm isfinite như sau: end isfinite(NaN) = 0 end isfinite(inf)=0 clear Tmin1 lamdamax Tmax0 isfinite(Số bất kỳ)=1) Tmin1=mang1 % Khu bo NaN trong mang Tmin1 lamdamax=mang1a A=isfinite(Tmin1) Tmax0=mang1b dem1=1; % Khu bo NaN trong mang Tmax1 for bien1=1:length(Tmin1) B=isnan(Tmax1) if (A(bien1)==1) dem2=1; mang1(dem1)=Tmin1(bien1); for bien2=1:length(Tmax1) mang1a(dem1)=lamdamax(bien1); if (B(bien2)==0) mang1b(dem1)=Tmax0(bien1); mang2(dem2)=Tmax1(bien2); dem1=dem1+1; mang2a(dem2)=lamdamin(bien2); end mang2b(dem2)=Tmin0(bien2); end dem2=dem2+1; clear Tmin1 lamdamax Tmax0 end Tmin1=mang1 end lamdamax=mang1a clear Tmax1 lamdamin Tmin0 VIC FAN
SWEET NOVEMBER 29 VERSION 2009 Tmax1=mang2 for u=1:3 lamdamin=mang2a Tmin(length(Tmin))=[]; Tmin0=mang2b end 4 mảng thu được lúc này là: Tmin lamdamax = % Thong nhat mang Tmax0 va mang Tmax1 410 426 443 462 484 510 537 569 604 thanh mang Tmax, sap xep cac phan tu 645 692 % theo thu tu tuong ung voi buoc song tang lamdamin = dan 416 433 452 472 496 523 552 667 719 for v=1:min(length(Tmax1),length(Tmax0)) Tmin0 = Tmax(2*v-1)=Tmax0(v); 62.7780 66.6750 71.0550 72.4700 71.7250 Tmax(2*v)=Tmax1(v); 74.3590 75.1270 75.7660 74.3490 end Tmax0 = for u=1:3 69.2940 74.8510 79.3090 81.1100 84.9150 Tmax(length(Tmax))=[]; 86.1450 86.3630 86.7180 86.9350 86.9110 end 87.3750 Tmax Đến đây thì tạm ổn, mặc dù “thiên hạ vẫn Kết quả: chưa thái bình”. Bước kế tiếp bạn cần quan tâm lamda = là thống nhất 2 mảng bước sóng lamdamin và 410 416 426 433 443 452 462 472 484 lamdamax thành một mảng lamda duy nhất; kết 496 510 523 537 552 569 Tmin = hợp 2 mảng Tmin0 và Tmin1 thành một mảng 60.1236 62.7780 65.0704 66.6750 68.9803 Tmin duy nhất; nối kết 2 mảng Tmax0 và 71.0550 71.7625 72.4700 72.0975 71.7250 Tmax1 thành một mảng Tmax duy nhất. Nếu 73.0908 74.3590 74.7298 75.1270 75.2215 bạn bỏ qua bước này thì xem như những yêu Tmax = cầu ở sau bạn hoàn toàn bế tắc. 69.2940 71.3779 74.8510 76.6866 79.3090 Ở đây, tác giả cung cấp cho bạn 2 cách 80.1621 81.1100 82.8395 84.9150 85.4827 thống nhất mảng, tùy bạn đọc chọn 1 mà dùng. 86.1450 86.2500 86.3630 86.5294 86.7180 Cách 1: Chỉ áp dụng cho bài này và bảng số Cách 2: liệu này, chưa chính xác lắm, các kết quả thu % Thong nhat cac mang được dựa trên cách 2 và cách 3 lamda=[lamdamin lamdamax]; % Thong nhat mang lamdamin va lamdamax Tmin=[Tmin0 Tmin1]; thanh mang lamda Tmax=[Tmax1 Tmax0]; for % Sap xep cac mang theo thu tu buoc song v=1:min(length(lamdamin),length(lamdamax)) tang dan lamda(2*v-1)=lamdamax(v); for v=1:length(lamda)-1 lamda(2*v)=lamdamin(v); for u=v+1:length(lamda) end if (lamda(v)>lamda(u)) for u=1:3 tam1=lamda(v); lamda(length(lamda))=[]; lamda(v)=lamda(u); end lamda(u)=tam1; lamda tam2=Tmin(v); % Thong nhat mang Tmin0 va mang Tmin1 Tmin(v)=Tmin(u); thanh mang Tmin, sap xep cac phan tu Tmin(u)=tam2; % theo thu tu tuong ung voi buoc song tang tam3=Tmax(v); dan Tmax(v)=Tmax(u); for v=1:min(length(Tmin1),length(Tmin0)) Tmax(u)=tam3; Tmin(2*v-1)=Tmin1(v); end Tmin(2*v)=Tmin0(v); end end end lamda VIC FAN
SWEET NOVEMBER 30 VERSION 2009 Tmin % Chiet suat cua de: Tmax S=1.52; Cách 3: % Xac dinh N tuong ung voi moi cap gia tri % Thong nhat cac mang, sap xep cac mang Tmax va Tmin theo thu tu buoc song tang dan % (Ghi chu: Cac gia tri cua TM va Tm duoc doi lamda_T=[lamdamin lamdamax;Tmin0 Tmin1; tu % ra so thap phan) Tmax1 Tmax0]; for o=1:length(lamda) lamda=[lamdamin lamdamax]; N(o)=2*S*(Tmax(o)/100-Tmin(o)/100) for x=1:length(lamda)-1 /((Tmax(o)/100)*(Tmin(o)/100))+(S^2+1)/2; for y=x+1:length(lamda) end if (lamda(x)>lamda(y)) N gan=lamda_T(:,x); % Xac dinh chiet suat n1 cua mang tuong ung lamda_T(:,x)=lamda_T(:,y); voi moi cap gia tri TM va Tm lamda_T(:,y)=gan; for p=1:length(lamda) lamda=lamda_T(1,:); n1(p)=sqrt(N(p)+(N(p)^2-S^2)^(1/2)); end end end n1 end Đây là kết quả bạn thu được: lamda=lamda_T(1,:) N = Tmin=lamda_T(2,:) 2.3243 2.2386 2.2657 2.2504 2.2291 Tmax=lamda_T(3,:) 2.1413 2.1434 2.1803 2.2917 2.3373 Đến đây thì kết cục “đại đoàn viên”. Bạn 2.2855 2.2188 2.2032 2.1884 2.1910 thử xuất ra các mảng lamda, lamdamin, 2.1893 2.1762 2.1784 2.2247 2.2531 lamdamax xem sao n1 = lamda = 2.0206 1.9703 1.9864 1.9774 1.9646 410 416 426 433 443 452 462 472 1.9103 1.9117 1.9348 2.0017 2.0280 484 496 510 523 537 552 569 604 1.9981 1.9584 1.9488 1.9398 1.9414 645 667 692 719 1.9404 1.9323 1.9336 1.9619 1.9789 Tmin = Chúng ta bắt tay làm tiếp câu 3 60.1236 62.7780 65.0704 66.6750 % Cau 3 68.9803 71.0550 71.7625 72.4700 % Xac dinh do day di cua mang tuong ung voi 72.0975 71.7250 73.0908 74.3590 tat ca cac cuc tri 74.7298 75.1270 75.2215 75.4159 for q=1:(length(lamda)-1) 75.6438 75.7660 75.0848 74.3490 tu(q)=lamda(q)*lamda(q+1); Tmax = mau(q)=2*(lamda(q+1)*n1(q)- 69.2940 71.3779 74.8510 76.6866 lamda(q)*n1(q+1)); 79.3090 80.1621 81.1100 82.8395 doday(q)=tu(q)/mau(q); 84.9150 85.4827 86.1450 86.2500 end doday % Lay ca gia tri am 86.3630 86.5294 86.7180 86.9350 Kết quả: 86.9110 87.1282 87.3750 87.0827 doday = Thế thì ổn rồi! Vậy là phần khó nhất của 1.0e+004 * bài 1 cũng xong xuôi. Các câu còn lại vô cùng 0.2605 0.6810 0.5195 0.3791 0.2400 dễ dàng vì bạn chỉ cần ghép vào công thức là 0.5646 1.2898 1.3675 1.0660 0.2924 xong. (Phần này chắc ai cũng làm giống ai, 0.2886 0.4334 0.4348 0.4891 0.2507 bạn đọc có thể thay đổi ký hiệu các biến hoặc 0.2307 0.5169 0.7831 0.6039 tham khảo thêm các bài của các anh chị đi % Cau 4: trước) % Cach 1: Dung ham % Cau 2 % Tinh do day d1 cua cac di vua tim d1=mean(doday) VIC FAN
SWEET NOVEMBER 31 VERSION 2009 % Cach 2: Dung lap trinh for r=2:length(bacm) % Tinh do day d1 cua cac di vua tim m(r)=m(dem)+0.5; tong=0; dem=dem+1; for r=1:length(lamda)-1 end tong=tong+doday(r); for k1=1:length(m)-1 end for k2=k1+1:length(m) d1=tong/(length(lamda)-1); if (m(k1)<m(k2)) d1 k3=m(k1); Kết quả: m(k1)=m(k2); d1 = m(k2)=k3; 5.6272e+003 end % Cau 5: end % Cach 1 end % Dung cac chiet suat n va d1 de xac dinh bac m m cua cac cuc tri Ket qua: % Doi voi cac dinh cuc tieu tren pho, bac m la m = so ban nguyen 40.5000 40.0000 39.5000 39.0000 38.5000 for s=1:2:length(lamda) 38.0000 37.5000 37.0000 36.5000 36.0000 bacm(s)=(2*n1(s)*d1)/lamda(s); 35.5000 35.0000 34.5000 34.0000 33.5000 bac(s)=round(bacm(s)); 33.0000 32.5000 32.0000 31.5000 31.0000 if (bacm(s)<bac(s)) % Cach 3 bacm(s)=bac(s)-0.5; % Dung cac chiet suat n va d1 de xac dinh bac else m cua cac cuc tri bacm(s)=bac(s)+0.5; % Doi voi cac dinh cuc tieu tren pho, bac m la end so ban nguyen end for s=1:length(lamda) % Doi voi cac dinh cuc dai tren pho, bac m la m(s)=(2*n1(s)*d1)/lamda(s); so nguyen end m(length(m))=round(m(length(m))) for t=2:2:length(lamda) for r=length(m)-1:-1:1 bacm(t)=(2*n1(t)*d1)/lamda(t); m(r)=m(r+1)+0.5; bacm(t)=round(bacm(t)); end end m bacm Ket qua: Ket qua: m = bacm = 47.0000 46.5000 46.0000 45.5000 45.0000 53.5000 51.0000 50.5000 50.0000 48.5000 44.5000 44.0000 43.5000 43.0000 42.5000 46.0000 44.5000 44.0000 44.5000 44.0000 42.0000 41.5000 41.0000 40.5000 40.0000 42.5000 41.0000 39.5000 38.0000 37.5000 23.0000 35.0000 32.5000 31.0000 30.5000 30.0000 % Cau 6: 28.5000 27.0000 25.5000 25.0000 23.5000 % Tinh cac do day di2 23.0000 for u=1:length(lamda) % Cach 2 di2(u)=m(u)*lamda(u)/(2*n1(u)); % Dung cac chiet suat n va d1 de xac dinh bac end m cua cac cuc tri di2 % Doi voi cac dinh cuc tieu tren pho, bac m la % Tinh do day trung binh d2 so ban nguyen % Cach 1 for s=1:length(lamda) tongd=0; bacm(s)=(2*n1(s)*d1)/lamda(s); for v=1:length(lamda) end tongd=tongd+di2(v); m=round(min(bacm)); end dem=1; d2=tongd/(length(lamda)); VIC FAN
SWEET NOVEMBER 32 VERSION 2009 d2 % Cach 2 d2=mean(di2) Ket qua: di2 = 1.0e+003 * 4.1089 4.2227 4.2355 4.2701 4.3407 4.4955 4.5313 4.5132 4.4128 4.4023 4.5307 4.6735 4.7532 4.8376 4.9093 5.1362 5.4244 5.5191 5.5552 5.6315 d2 = 4.7252e+003 % Cau 7 % Tim lai chiet suat n2 for w=1:length(lamda) n2(w)=(m(w)*lamda(w))/(2*d2); end n2 Ket qua n2 = 1.7571 1.7608 1.7806 1.7869 1.8047 1.8175 1.8333 1.8480 1.8693 1.8895 1.9158 1.9370 1.9604 1.9860 2.0170 2.1091 2.2182 2.2585 2.3066 2.3585 Xong xuôi, bạn đọc ghi lại các kết quả cần thiết vào bảng số liệu: BÀI LẬP TRÌNH 2 clc Trước khi vào bài lập trình 2, mời các bạn clear all làm quen với một số vấn đề. fprintf('Nghiem cua cac phuong trinh la:') 1) Dùng những hàm có sẵn trong Matlab, hãy % Cau a giải các phương trình sau: nghiem1=solve('2*x-3=0') % Cau b a> 2x – 3 = 0 nghiem2=solve('89+12.5*x=0') b> 89 + 12.5x = 0 % Cau c 2 c> x 19.77 0 nghiem3=solve('2/3*x-19.77=0') 3 % Cau d d> x2 + 2x + 1 = 0 nghiem4=solve('x^2+2*x+1=0') e> -x2 – 2x + 3 = 0 % Cau e f> 3x2 + 3 = 0 nghiem5=solve('-x^2-2*x+3=0') g> 2x3 + 5x2 – 3x – 4 = 0 % Cau f nghiem6=solve('3*x^2+3=0') Hướng dẫn: Để giải phương trình hay hệ % Cau g phương trình trong Matlab, bạn có thể sử dụng nghiem7=solve('2*x^3+5*x^2-3*x-4=0') hàm solve với cấu trúc tổng quát như sau: solve(‘phương trình 1’, ‘phương trình 2’, , ‘phương trình N’) solve(‘phương trình 1’, ‘phương trình 2’, , ‘phương trình N’, ‘bien1’, ‘bien2’, ,‘bien N’) Trong M-file, bạn gõ các dòng chương trình sau: VIC FAN
SWEET NOVEMBER 33 VERSION 2009 Tiếp tục dùng hàm solve, chúng ta sẽ tiến hành giải các hệ phương trình: clc clear all fprintf('Nghiem cua cac he phuong trinh la:'); % Cau a [nghiemx1,nghiemy1]=solve('2*x-3*y-1=0','- 5*x+4*y+6=0','x','y') % Cau b [nghiemx2,nghiemy2]=solve('4*(x^2-2*x)- 3*(y^2-y)=-10','5*(x^2-2*x)-2*(y^2-y)=-1','x','y') % Cau c [nghiemx3,nghiemy3]=solve('1/x+8/y=17','7/x- 3/y=1','x','y') 2) Dùng những hàm có sẵn trong Matlab giải các hệ phương trình sau: 2x 3y 1 0 a> 5x 4y 6 0 4 x2 2x 3 y2 y 10 b> 3) Quang ma trận (Matrix Optics) 2 2 5 x 2x 2 y y 1 Quang ma trận là kỹ thuật tính toán áp dụng 1 8 ma trận để mô tả sự lan truyền của tia sáng qua 17 các thành phần quang học (như không gian tự x y c> do, mặt cầu, mặt gương, thấu kính, gương 7 3 1 cầu, ), qua đó thiết lập mối liên hệ giữa tia x y sáng tới và tia sáng ló. VIC FAN
SWEET NOVEMBER 34 VERSION 2009 Giả sử một tia sáng đi vào một hệ quang học 1 0 A B tại vị trí y với góc θ và đi ra hệ quang học tại n 1 1 M = 0 1 vị trí y2 với góc θ 2. Cặp giá trị (y2,θ2) có thể C D n2 biểu diễn theo cặp giá trị (y ,θ ) thông qua hệ 2 1 1  Khúc xạ tại mặt phân cách là mặt cầu bán phương trình: kính R y2 Ay1 B1 2 Cy1 D1 Với A, B, C, D là các số thực. Hệ 2 phương trình trên có thể được viết lại dưới dạng ma trận: y2 A B y1 2 C D 1 A B Trong đó ma trận M = được gọi là C D ma trận truyền tia (The ray-transfer matrix), A B 1 0 n n được xem là đặc trưng cho hệ quang học. M = 2 1 C D 1 Ma trận truyền tia của một số thành phần R quang học đơn giản:  Ma trận phản xạ  Không gian tự do  Phản xạ tại gương phẳng n d A B 1 M = n A B 1 0 C D 0 1 M = C D 0 1  Ma trận khúc xạ  Phản xạ tại gương cầu  Khúc xạ tại mặt phân cách là mặt phẳng VIC FAN
SWEET NOVEMBER 35 VERSION 2009 A B 1 0 1 d / n M = 2 Môi trường chiết suất n: Mkk = C D 1 0 1 R Mặt phân cách cầu:  Thấu kính mỏng tiêu cự f A B 1 0 Mmc = n2 n1 C D 1 R Hướng dẫn: Hệ quang học của chúng ta gồm 5 thành phần quang học được xếp theo thứ tự: Không khí M1 → Mặt cầu phân cách bán kính r1 M2 → Môi trường thấu kính M 3 → Mặt cầu phân cách bán kính r2 M4 → Không khí M5 A B 1 0 a) Trong M-file, bạn gõ các dòng lệnh sau: 1 M = clc C D 1 f clear all Ma trận truyền tia qua một hệ quang học % Nhap cac du lieu gồm nhiều thành phần quang học: n1=input('Nhap vao chiet suat khong khi:'); n2=input('Nhap vao chiet suat thau kinh:'); Tia tới → M → M → → M → M → Tia ló 1 2 n-1 n BC=input('Nhap vao be day thau kinh:'); A B r1=input('Nhap vao ban kinh thu nhat cua thau M = Mn.Mn-1.Mn-2 M3.M2.M1 = C D kinh:'); Điều kiện để có ảnh: B = 0 r2=input('Nhap vao ban kinh thu hai cua thau kinh:'); % Cau a 4) Một số bài tập mô phỏng hệ quang ma trận syms CD h2 % Khai bao bien su dung la CD va Bài tập 1: Cho hệ quang học như hình vẽ: h2 AB=input('Nhap vao khoang cach giua vat va n1 thau kinh:'); h1=input('Nhap vao chieu cao vat:'); % Cac ma tran dac trung cho cac thanh phan quang hoc n2 M1=[1 AB/n1;0 1]; M2=[1 0;-(n2-n1)/r1 1]; M3=[1 BC/n2;0 1]; M4=[1 0;(n1-n2)/r2 1]; Thấu kính có bề dày BC = 0.5cm, chiết suất M5=[1 CD/n1;0 1]; M=M5*M4*M3*M2*M1 n2 = 1.52 được tạo bởi hai mặt cầu bán kính lần lượt là r = 10cm, r = 20cm được đặt trong A=M(1,1) 1 2 B=M(1,2) không khí có chiết suất n1 = 1. C=M(2,1) a) Vật có chiều cao h 1 = 1cm, cách thấu kính D=M(2,2) một khoảng AB = 5cm. Hãy xác định vị trí CD CD=double(solve(B)) % solve(B): Giai tim CD; và chiều cao h2 của ảnh. double( ): Chuyen ket qua sang so thap phan b) Giả sử biết ảnh cách thấu kính một khoảng % Ta co h2=A*h1+B*theta1 (xem trang 31); vi CD = 10cm, chiều cao h2 = 2cm. Hãy xác định B=0 (dieu kien de co anh) nen h2=A*h1 lại vị trí AB và chiều cao h1 của vật. h2=subs(A*h1) % A*h1= Bieu thuc chua CD; Cho ma trận truyền tia của các thành phần subs( ): The CD vua tim duoc vao bieu thuc de như sau: tinh h2 Các kết quả xuất ra: VIC FAN
SWEET NOVEMBER 36 VERSION 2009 A = h1=subs(A*h2) % A*h2= Bieu thuc chua AB; 747/760-29471/380000*CD subs( ): The AB vua tim duoc vao bieu thuc de tinh B = h1 797/152+45879/76000*CD Các kết quả xuất ra: C = A = -29471/380000 1507/1520-29471/380000*AB D = B = 45879/76000 1557/152+7879/38000*AB CD = C = -8.6859 h2 = -29471/380000 1.6565 D = (Nhận xét: Dấu trừ ở CD thể hiện ảnh nằm 7879/38000 AB = cùng phía và cùng chiều với vật với độ cao h 2 = 1.6565cm) -49.4035 b) Tương tự như câu a nhưng lúc này các ma h1 = trận có sự thay đổi. Bạn đọc tham khảo đoạn 9.6459 chương trình bên dưới: clc BÀI LẬP TRÌNH 3 clear all XÁC ĐỊNH PHỔ TRUYỀN QUA CỦA % Nhap cac du lieu MÀNG DỰA VÀO BƯỚC SÓNG, ĐỘ DÀY n1=input('Nhap vao chiet suat khong khi:'); VÀ PHƯƠNG TRÌNH CHIẾT SUẤT n2=input('Nhap vao chiet suat thau kinh:'); Các yêu cầu trong bài lập trình này: BC=input('Nhap vao be day thau kinh:'); 1) Lập trình chọn vùng bước sóng hoạt động, r1=input('Nhap vao ban kinh thu nhat cua thau một trong 3 vùng: vùng truyền suốt (0 – kinh:'); 0.4μm), vùng hấp thụ yếu và trung bình (0.4 – r2=input('Nhap vao ban kinh thu hai cua thau 0.7μm), vùng hấp thụ mạnh (0.7 – 1μm). kinh:'); % Cau b 2) Lập trình chọn bước sóng nhỏ nhất và lớn syms AB h1 % Khai bao bien su dung la AB va h1 nhất tương ứng với các vùng đã chọn ở trên. CD=input('Nhap vao khoang cach giua anh va 3) Tạo mảng bước sóng. Chọn bước nhảy của thau kinh:'); bước sóng là 2nm = 0.002μm. h2=input('Nhap vao chieu cao anh:'); 4) Viết chương trình chọn chất cần mô phỏng, % Cac ma tran dac trung cho cac thanh phan một trong số 12 chất sau: PbTe, Ge, CdTe, quang hoc ZnSe, ZnS, Ta2O5, YbF3, YF3, SiO2, MgF2, M1=[1 CD/n1;0 1]; CaF2, BaF2. M2=[1 0;-(n2-n1)/r2 1]; 5) Từ phương trình chiết suất của 12 chất, hãy M3=[1 BC/n2;0 1]; tính chiết suất của từng chất thay đổi theo bước M4=[1 0;(n1-n2)/r1 1]; sóng. Cho biết phương trình chiết suất cụ thể M5=[1 AB/n1;0 1]; M=M5*M4*M3*M2*M1 của từng chất là: A=M(1,1)  PbTe (Lead Telluride) B=M(1,2) 8.0749 107 n 5.0781 C=M(2,1) 4.  6715 2 6.0882 107 D=M(2,2) AB=double(solve(B)) % solve(B): Giai tim AB;  Ge (Germanium) double( ): Chuyen ket qua sang so thap phan 6.0528 107 n 4.0517 % Ta co h1=A*h2+B*theta2; vi B=0 (dieu kien de  8131.5 2 7.9566 106 co anh) nen h1=A*h2  CdTe (Cadmium Telluride) VIC FAN
SWEET NOVEMBER 37 VERSION 2009 9.2789 107 A = 16n2S n 2.0681 3 2  8157.6 2 1.2387 108 B = (n+1) (n+S ) C = 2(n2-1)(n2-S2)  ZnSe (Zinc Selenide) D = (n-1)3(n-S2) 3.4469 108 n 2.0845 φ = 4πnd/λ  6715 2 2.7301 108 x = exp(-αd)  ZnS (Zinc Sulphide) 8) Vẽ đồ thị độ truyền qua T thay đổi theo bước . 7.0775 108 sóng của chất Ta2O5 n 1.7951 2 8 Hướng dẫn:  35622 8.3564 10 Đây là bài toán ngược của bài lập trình 1  Ta2O5 (Tantalum Pentoxide) nhưng lại ở mức độ đơn giản hơn và dễ lập 6.4378 107 trình hơn. Bốn câu đầu bạn có thể dễ dàng đối n 1.9079  27923 2 2.2668 109 phó mà không chút quá khó. Bạn có thể tham khảo vài đoạn chương trình dưới đây.  YbF3 (Ytterbium Fluoride) 9 CÂU 1 2.7005 10 Cách 1: n 20.09287 2  15717 1.9486 109 clc  YF3 (Yttrium Fluoride) clear all 4.8652 108 % Cau 1: Lap trinh chon 1 trong 3 vung hoat dong n 1.0527 2 9  14675 1.1942 10 fprintf('Chon vung hoat dong. Nhap so 1 ung  SiO2 (Silicon Dioxide) voi vung truyen suot, so 2 voi vung hap thu yeu 4.3456 109 - trung binh, so 3 voi vung hap thu manh. ') key1=input('Nhap:'); n 0.76991 2 9  40506 4.4229 10 if (key1~=1)&(key1~=2)&(key1~=3)  MgF2 (Magnesium Fluoride) error('So lieu khong hop le. Vui long nhap 1.9465 108 lai.') n 1.371 2 9 end  132800 2.9588 10 Cách 2: Chỉ cần thay fprintf thành disp, bạn đã  CaF2 (Calcium Fluoride) có thêm 1 cách mới. 1.1747 109 clc n 1.1199  41379 2 4.8958 109 clear all % Cau 1: Lap trinh chon 1 trong 3 vung hoat  BaF2 (Barium Fluoride) dong 1.50799 109 disp('Chon vung hoat dong. Nhap so 1 ung voi n 0.86148  30788 2 2.3347 109 vung truyen suot, so 2 voi vung hap thu yeu - trung binh, so 3 voi vung hap thu manh.') 1 Trong đó:  key1=input('Nhap:');  if (key1~=1)&(key1~=2)&(key1~=3) 6) Viết chương trình chọn một trong 3 độ dày error('So lieu khong hop le. Vui long nhap màng: 0.5μm, 5μm và 15μm. lai.') end 7) Chọn Ta2O5 làm chất khảo sát. Tính độ truyền qua theo bước sóng dựa vào công thức Ở đây, chúng ta thấy có sự xuất hiện của (3.8): hàm error. Hàm này có cái độc đáo là nếu CT (3.8): người chạy chương trình phạm lỗi nhập sai, nó Ax sẽ báo lỗi và vứt bỏ nhiệm vụ, không làm tiếp T nữa. Hay chưa? B Cxcos Dx2 Trong đó: VIC FAN
SWEET NOVEMBER 38 VERSION 2009 Bạn cũng có thể sử dụng vòng lặp, bắt người Cách 2: Dùng vòng lặp while chạy nhập hoài nhập hoài cho đến khi số liệu % Cau 2: Chon buoc song nho nhat, lon nhat hợp lệ thì thôi. Cách 3 dưới đây là một thí dụ, tuong ung voi tung vung cho thấy hàm while được sử dụng hiệu quả thế lamdamin=input('Nhap vao buoc song nho nào. nhat:') Cách 3: lamdamax=input('Nhap vao buoc song lon nhat:') clc while (lamdamax 0.4)) trung binh, so 3 voi vung hap thu manh.') disp('Vung truyen suot') key1=input('Nhap:'); lamdamin=input('Nhap lai buoc song nho while (key1~=1)&(key1~=2)&(key1~=3) nhat:') key1=input('So lieu khong hop le. Vui long lamdamax=input('Nhap lai buoc song lon nhap lai. Nhap:'); nhat:') end end CÂU 2 while “Phóng lao thì phải theo lao”, vì có 3 vùng (key1==2)&((lamdamin 0.7)) phổ hoạt động nên giờ đây bạn cũng có 3 mảng disp('Vung hap thu yeu - trung binh') bước sóng khác nhau. Cách lập trình khéo léo lamdamin=input('Nhap lai buoc song nho bây giờ là điều cần thiết. Dưới đây là vài cách nhat:') tham khảo. lamdamax=input('Nhap lai buoc song lon Cách 1: Dùng lenh dieu kien if nhat:') % Cau 2: Chon buoc song nho nhat, lon nhat end tuong ung voi tung vung while lamdamin=input('Nhap vao buoc song nho (key1==3)&((lamdamin 1)) nhat:') disp('Vung truyen suot') lamdamax=input('Nhap vao buoc song lon lamdamin=input('Nhap lai buoc song nho nhat:') nhat:') if (lamdamax =0.35)&(lamdamax =0.4)&(lamdamax =0.7)&(lamdamax<=1) if (lamdamax<=lamdamin) fprintf('Vung hap thu manh') error('Vui long nhap lai, lamdamax phai lon else hon lamdamin') error('Khong xac dinh duoc vung pho. Vui end long kiem tra lai lamdamin va lamdamax') end VIC FAN
SWEET NOVEMBER 39 VERSION 2009 while Cách 1: (key1==1)&((lamdamin 0.4)) % Cau 3: Tao mang buoc song, buoc nhay disp('Vung truyen suot') 0.002 micromet lamdamin=input('Nhap lai buoc song nho dem=1; nhat:') buocnhay=0.002; lamdamax=input('Nhap lai buoc song lon for lamda=lamdamin:buocnhay:lamdamax nhat:') buocsong(dem)=lamda; end xicma(dem)=1/lamda; while dem=dem+1; (key1==2)&((lamdamin 0.7)) end disp('Vung hap thu yeu - trung binh') buocsong lamdamin=input('Nhap lai buoc song nho xicma nhat:') Cách 2: lamdamax=input('Nhap lai buoc song lon % Cau 3: Tao mang buoc song, buoc nhay nhat:') 0.002 micromet end dem=1; while buocnhay=0.002; (key1==3)&((lamdamin 1)) lamda=lamdamin; disp('Vung truyen suot') while (lamda =0)&(lamdamax =0.4)&(lamdamax =0.7)&(lamdamax 12) VIC FAN
SWEET NOVEMBER 40 VERSION 2009 error('So lieu nhap vao khong hop le. Vui d=[6715 8131.5 8157.6 6715 35622 27923 long kiem tra lai.') 15717 14675 40506 132800 41379 30788]; end f=[6.0882*10^7 7.9566*10^6 1.2387*10^8 Cách 2: 2.7301*10^8 8.3564*10^8 2.2668*10^9 % Cau 4: Chon chat can mo phong trong 12 1.9486*10^9 1.1942*10^9 4.4229*10^9 chat 2.9588*10^9 4.8958*10^9 2.3347*10^9]; disp('Chon chat can mo phong theo ma so:') % Tinh chiet suat theo su thay doi cua buoc disp('So 1: PbTe, So 2: Ge, So 3: CdTe, So 4: song theo cac phuong trinh chiet suat ZnSe, So 5: ZnS, So 6: Ta2O5') % Xuong dong for g=1:length(xicma) cho de nhin mau(g)=c(key2)*(xicma(g)- disp('So 7: YbF3, So 8: YF3, So 9: SiO2, So d(key2))^2+f(key2); 10: MgF2, So 11: CaF2, So 12: BaF2') chietsuat(g)=a(key2)+b(key2)/mau(g); key2=input('Nhap:') end while (key2 12) chietsuat key2=input('Chua dung. Vui long nhap lai:') Cách 2: end % Cau 5: Tim phuong trinh chiet suat tong quat CÂU 5 % Bieu thuc tinh chiet suat cua 12 chat dau Dễ dàng nhận thấy rằng phương trình chiết tien co dang tong quat suất của 12 chất đề bài cho đều có nét hau hau % n=a+b/(c*(xicma-d)^2+f) tương đồng, “tâm đầu ý hợp”. Thế thì tại sao % a,b,c,d,f la cac hang so dac trung cho moi bạn không làm một phương trình tổng quát đại chat % Ta se lap cac mang du lieu a,b,c,d,f nhu b sau: khái như là n a 2 để dễ dàng c  d f a=[5.0781 4.0517 2.0681 2.0845 1.7951 “thao túng” cả 12 em. Bạn sẽ thành lập trực 1.9079 0.90287 1.0527 0.76991 1.371 1.1199 tiếp 5 mảng dữ liệu a, b, c, d, f. Tất nhiên nếu là 0.86148]; 12 chất thì sẽ là mảng 12 phần tử. Bạn khéo léo b=[8.0749*10^7 6.0528*10^7 9.2789*10^7 và từ tốn lập trình cẩn thận và thiệt hay phần 3.4469*10^8 7.0775*10^8 6.4378*10^8 2.7005*10^9 4.8652*10^8 4.3456*10^9 mảng chiết suất theo mảng sigma. Đến đây thì 1.9465*10^8 1.1747*10^9 1.50799*10^9]; không cần biết bên trên người chạy chọn chất c=[4 1 1 1 1 1 1 1 1 1 1 1]; gì, bạn luôn có sẵn một mảng chiết suất theo d=[6715 8131.5 8157.6 6715 35622 27923 bước sóng (thật ra là theo sigma) tương ứng. 15717 14675 40506 132800 41379 30788]; Cách 1: f=[6.0882*10^7 7.9566*10^6 1.2387*10^8 % Cau 5: Tim phuong trinh chiet suat tong quat 2.7301*10^8 8.3564*10^8 2.2668*10^9 % Bieu thuc tinh chiet suat cua 12 chat dau 1.9486*10^9 1.1942*10^9 4.4229*10^9 tien co dang tong quat 2.9588*10^9 4.8958*10^9 2.3347*10^9]; % n=a+b/(c*(xicma-d)^2+f) % Tinh chiet suat theo su thay doi cua buoc % a,b,c,d,f la cac hang so dac trung cho moi song theo cac phuong trinh chiet suat chat for g=1:length(xicma) % Ta se lap cac mang du lieu a,b,c,d,f nhu chietsuat(g)=a(key2)+b(key2)/( sau: c(key2)*(xicma(g)-d(key2))^2+f(key2)); a=[5.0781 4.0517 2.0681 2.0845 1.7951 end 1.9079 0.90287 1.0527 0.76991 1.371 1.1199 chietsuat 0.86148]; Đây là kết quả thu được: b=[8.0749*10^7 6.0528*10^7 9.2789*10^7 chietsuat = 3.4469*10^8 7.0775*10^8 6.4378*10^8 2.1192 2.1192 2.1192 2.1192 2.1192 2.7005*10^9 4.8652*10^8 4.3456*10^9 2.1192 2.1192 2.1192 2.1192 2.1192 1.9465*10^8 1.1747*10^9 1.50799*10^9]; 2.1192 2.1192 2.1192 2.1192 2.1192 c=[4 1 1 1 1 1 1 1 1 1 1 1]; 2.1192 2.1192 2.1192 2.1192 2.1192 VIC FAN
SWEET NOVEMBER 41 VERSION 2009 2.1192 2.1192 2.1192 2.1192 2.1192 doday=input('Nhap lai do day mang 2.1192 2.1192 2.1192 2.1192 2.1192 (micromet), chon 1 trong 3: 0.5, 5, 15:') 2.1192 2.1192 2.1192 2.1192 2.1192 end 2.1192 2.1192 2.1192 2.1192 2.1192 CÂU 7 2.1192 2.1192 2.1192 2.1192 2.1192 Cách 1: 2.1192 2.1192 2.1192 2.1192 2.1192 % Cau 7: Tinh do truyen qua theo buoc song 2.1192 2.1192 2.1192 2.1192 2.1192 % Chiet suat cua de: S=1.52; 2.1192 2.1192 2.1192 2.1192 2.1192 for l=1:size(buocsong,2) 2.1192 2.1192 2.1192 2.1192 2.1192 A(l)=16*chietsuat(l)^2*S; 2.1192 2.1192 2.1192 2.1192 2.1192 B(l)=(chietsuat(l)+1)^3*(chietsuat(l)+S^2); 2.1192 2.1192 2.1192 2.1192 2.1192 C(l)=2*(chietsuat(l)^2-1)*(chietsuat(l)^2- 2.1192 2.1192 2.1192 2.1192 2.1192 S^2); 2.1192 2.1192 2.1192 2.1192 2.1192 D(l)=(chietsuat(l)-1)^3*(chietsuat(l)-S^2); 2.1192 2.1192 2.1192 2.1192 2.1192 phi(l)=4*pi*chietsuat(l)*doday/buocsong(l); 2.1192 2.1192 2.1192 2.1192 2.1192 alpha(l)=4*pi*chietsuat(l)*doday/buocsong(l); 2.1192 2.1192 2.1192 2.1192 2.1192 x(l)=exp(-alpha(l)*doday*10^-6); 2.1192 2.1192 2.1192 2.1192 2.1192 T(l)=A(l)*x(l)/(B(l)- 2.1192 2.1192 2.1192 2.1192 2.1192 C(l)*x(l)*cos(phi(l))+D(l)*x(l)^2); end 2.1192 2.1192 2.1192 2.1192 2.1192 T 2.1192 2.1192 2.1192 2.1192 2.1192 Cách 2: 2.1192 2.1192 2.1192 2.1192 2.1192 % Cau 7: Tinh do truyen qua theo buoc song 2.1192 2.1192 2.1192 2.1192 2.1192 S=1.52; 2.1192 2.1192 2.1192 2.1192 2.1192 for l=1:size(buocsong,2) 2.1192 2.1192 2.1192 2.1192 2.1192 A=16*chietsuat(l)^2*S; 2.1192 2.1192 2.1192 2.1192 2.1192 B=(chietsuat(l)+1)^3*(chietsuat(l)+S^2); 2.1192 2.1192 2.1192 2.1192 2.1192 C=2*(chietsuat(l)^2-1)*(chietsuat(l)^2-S^2); Trời, thế này là thế nào??? Bạn đừng vội D=(chietsuat(l)-1)^3*(chietsuat(l)-S^2); hoang mang hay nghĩ rằng mình sai. Đó là điều phi=4*pi*chietsuat(l)*doday/buocsong(l); xảy ra khi “lấy trứng chọi đá”, ý là một số quá alpha=4*pi*chietsuat(l)*doday/buocsong(l); nhỏ đối đầu với một số quá lớn, kết quả hiển thị x=exp(-alpha*doday*10^-6); T(l)=A*x/(B-C*x*cos(phi)+D*x^2); trong chừng mực không cho thấy sự khác biệt. end Chúng ta cứ tiếp tục “lên đường” thì sẽ rõ thực T hư. Kết quả thu được với vùng hấp thụ yếu và CÂU 6 trung bình, lamdamin = 0.4, lamdamax = 0.7 Cách 1: chất mô phỏng là Ta2O5 và độ dày là 5: % Cau 6: Chon do day cua mang T = doday=input('Nhap vao do day mang 0.9157 0.7943 0.7329 0.8436 0.9090 (micromet), chon 1 trong 3: 0.5, 5, 15:') if (doday~=0.5)&(doday~=5)&(doday~=15) 0.7795 0.7357 0.8487 0.9093 0.7849 error('So lieu nhap vao chua dung. Vui long 0.7327 0.8312 0.9159 0.8112 0.7304 nhap lai!') 0.7936 0.9095 end 0.8598 0.7464 0.7494 0.8630 0.9102 Cách 2: 0.8017 0.7305 0.7833 0.8974 0.8887 % Cau 6: Chon do day cua mang 0.7750 0.7311 0.8008 0.9058 0.8816 doday=input('Nhap vao do day mang 0.7722 0.7310 (micromet), chon 1 trong 3: 0.5, 5, 15:') while (doday~=0.5)&(doday~=5)&(doday~=15) VIC FAN
SWEET NOVEMBER 42 VERSION 2009 0.7946 0.8990 0.8945 0.7904 0.7309 0.7685 0.8707 0.9144 0.8340 0.7463 0.7380 0.8133 0.9047 0.8948 0.7996 0.7346 0.7495 0.8332 0.9117 0.8872 0.7945 0.7344 0.7474 0.8247 0.9059 0.8997 0.8156 0.7443 0.7350 0.7911 0.8780 0.9165 0.8629 0.7785 0.7324 0.7475 0.8160 0.8951 0.9127 0.8506 0.7719 0.7318 0.7473 0.8113 0.8884 0.9161 0.8675 0.7901 0.7385 0.7352 0.7804 0.8547 0.9115 0.9029 0.8380 0.7687 0.7325 0.7414 0.7919 0.8634 0.9132 0.9025 0.8413 0.7744 0.7352 0.7359 doday = 5micromet (xấu quá!) 0.7756 0.8409 0.9006 0.9152 0.8743 0.8083 0.7539 0.7307 0.7441 0.7900 0.8539 0.9058 0.9139 0.8731 0.8107 0.7578 0.7319 0.7384 0.7753 0.8329 0.8898 0.9170 0.8983 0.8459 0.7878 0.7459 0.7305 0.7438 0.7831 0.8386 0.8914 0.9169 0.9016 0.8546 0.7991 0.7548 0.7324 0.7353 0.7627 0.8093 0.8631 0.9051 0.9166 0.8922 0.8444 0.7928 0.7528 0.7324 0.7344 0.7580 0.7994 CÂU 8 (Quá dễ và đơn giản nên chỉ có 1 cách) % Ve do thi do truyen qua thay doi theo buoc doday = 15micromet (nguyên đám rừng luôn!!) song plot(buocsong,T) xlabel('Buoc song (micromet)') ylabel('Do truyen qua T(%)') title('Do thi do truyen qua thay doi theo buoc song') grid on Và sau đây là các kết quả bạn cần xuất ra: doday = 0.5micromet BÀI LẬP TRÌNH 4 Các yêu cầu trong bài lập trình này: VIC FAN
SWEET NOVEMBER 43 VERSION 2009 1) Lập trình chọn vùng bước sóng khảo sát Cách 1 (lamdamin, lamdamax, bước nhảy) clc 2) Lập trình chọn ra 2 chất bất kỳ trong số 12 clear all % 1 chất: PbTe, Ge, CdTe, ZnSe, ZnS, Ta2O5, buocsongmin=input(‘Nhap vao buocsongmin:’); YbF3, YF3, SiO2, MgF2, CaF2, BaF2. buocsongmax=input(‘Nhap vao buocsongmax:’); In kết quả với 2 chất là Ta2O5 và MgF2 buocnhay=input(“Nhap vao buoc nhay:’); 3) So sánh chiết suất của 2 chất vừa chọn ở câu buocsong=[buocsongmin:buocnhay:buocsongmax]; 2 tại phần tử đầu tiên của mảng bước sóng. Cách 2 Chất nào có chiết suất lớn hơn sẽ là chất chiết clc clear all suất cao, chất nào có chiết suất nhỏ hơn sẽ là % 1 chất chiết suất thấp. buocsongmin=input('Nhap vao buocsongmin:'); 4) Lập trình chọn chất đầu tiên là chất chiết buocsongmax=input('Nhap vao buocsongmax:'); suất cao hay chất chiết suất thấp. Chọn số lớp dem=1; màng cần mô phỏng. Sắp xếp xen kẽ các lớp buocnhay=input('Nhap buoc nhay:') for lamda=buocsongmin:buocnhay:buocsongmax màng theo thứ tự chiết suất cao – chiết suất buocsong(dem)=lamda; thấp – chiết suất cao – hay chiết suất thấp – xicma(dem)=1/lamda; chiết suất cao – chiết suất thấp – dem=dem+1; In kết quả với số lớp là 7. end 5) Nhập độ dày các lớp. buocsong In kết quả với độ dày (nm) theo thứ tự từ xa xicma Cách 3 đế đến gần đế là 245, 338, 81, 25, 30, 26, 146. clc 6) Tính độ phản xạ R của màng theo bước sóng clear all dựa vào công thức: % 1 B C 2 buocsongmin=input('Nhap vao buocsongmin:'); R  buocsongmax=input('Nhap vao buocsongmax:'); B C buocnhay=input('Nhap vao buoc nhay:'); Trong đó: lamda=buocsongmin; for bien=1:(buocsongmax-buocsongmin)/buocnhay+1 q isin j buocsong(bien)=lamda; B cos j 1 xicma(bien)=1/lamda; n j C  n lamda=lamda+buocnhay; j 1 in sin cos s j j j end 2 n d buocsong Với:  j j xicma m  Nói thêm: Độ truyền qua T của màng được CAU 2 tính theo công thức: Cách 1: disp('Chon ra 2 chat can mo phong theo ma so:') 4ns T  2 disp('So 1: PbTe So 2: Ge') B C disp('So 3: CdTe So 4: ZnSe') 7) Vẽ đồ thị độ phản xạ thay đổi theo bước disp('So 5: ZnS So 6: Ta2O5') sóng. disp('So 7: YbF3 So 8: YF3') disp('So 9: SiO2 So 10: MgF2') Hướng dẫn: disp('So 11: CaF2 So 12: BaF2') Đến bài này thì các bạn đã chính thức vào key1=input('Nhap chat thu nhat:'); “hang cọp”, tức loạt bài mô phỏng cực khó và if (key1 12) phức tạp nhất. Dưới đây là một số gợi ý, bạn error('So lieu nhap vao khong hop le. Vui long đọc nên đổi biến và tham khảo thêm bài khóa kiem tra lai.') end 02 và 03! key2=input('Nhap chat thu hai:'); CAU 1 if (key2 12) VIC FAN
SWEET NOVEMBER 44 VERSION 2009 error('So lieu nhap vao khong hop le. Vui long a=[5.0781 4.0517 2.0681 2.0845 1.7951 1.9079 kiem tra lai.') 0.90287 1.0527 0.76991 1.371 1.1199 0.86148]; end b=[8.0749*10^7 6.0528*10^7 9.2789*10^7 Cách 2: 3.4469*10^8 7.0775*10^8 6.4378*10^8 % 2 2.7005*10^9 4.8652*10^8 4.3456*10^9 disp('Ban hay chon ra 2 chat bat ky trong danh sach 1.9465*10^8 1.1747*10^9 1.50799*10^9]; 12 chat sau bang cach nhap vao ma so cua chat:') c=[4 1 1 1 1 1 1 1 1 1 1 1]; disp('PbTe: 1') d=[6715 8131.5 8157.6 6715 35622 27923 15717 disp('Ge: 2') 14675 40506 132800 41379 30788]; disp('CdTe: 3') f=[6.0882*10^7 7.9566*10^6 1.2387*10^8 disp('ZnSe: 4') 2.7301*10^8 8.3564*10^8 2.2668*10^9 disp('ZnS: 5') 1.9486*10^9 1.1942*10^9 4.4229*10^9 disp('Ta2O5:6') 2.9588*10^9 4.8958*10^9 2.3347*10^9]; disp('YbF3: 7') disp('Chiet suat cua chat thu nhat la:') disp('YF3: 8') for g1=1:length(xicma) disp('SiO2: 9') disp('MgF2:10') chietsuat1(g1)=a(key1)+b(key1)/(c(key1)*(xicma(g1 disp('CaF2:11') )-d(key1))^2+f(key1)); disp('BaF2:12') end key1= input('Nhap chat thu nhat:'); chietsuat1 if (key1 12) disp('Chiet suat cua chat thu hai la:') error('So lieu nhap vao khong hop le. Vui long for g2=1:length(xicma) kiem tra lai.') end chietsuat2(g2)=a(key2)+b(key2)/(c(key2)*(xicma(g2 key2=input('Nhap chat thu hai:'); )-d(key2))^2+f(key2)); if (key2 12) end error('So lieu nhap vao khong hop le. Vui long chietsuat2 kiem tra lai.') disp('Chat co chiet suat cao la chat thu:') end lon=max(chietsuat1(1),chietsuat2(1)); Cách 3: if (chietsuat1(1)==lon) % 2 key1 disp('Danh sach cac chat co the mo phong:') else fprintf('1-PbTe ') key2 fprintf('2-Ge ') end fprintf('3-CdTe ' ) disp('Chat co chiet suat thap la chat thu:') fprintf('4-ZnSe ') nho=min(chietsuat1(1),chietsuat2(1)); fprintf('5-ZnS ') if (chietsuat1(1)==nho) fprintf('6-Ta2O5 ') key1 fprintf('7-YbF3 ') else fprintf('8-YF3 ') key2 fprintf('9-SiO2 ') end fprintf('10-MgF2 ') fprintf('11-CaF2 ') CAU 4 fprintf('12-BaF2 ') % 4 disp('Nhap vao 2 chat') solop=input('Nhap vao so lop mang ban muon mo key1=input('Chat thu nhat:'); phong:') key2=input('Chat thu hai:'); disp('Cho biet lop mang dau tien la lop cua chat co while (key2==key1) chiet suat cao hay chiet suat thap.') key2=input(‘Chat thu hai trung voi chat thu nhat. disp(' Neu la chat co chiet suat cao thi nhap vao so Vui long nhap lai’) 1.') end disp(' Neu la chat co chiet suat thap thi nhap vao so 0.') CAU 3 tam=input(' Nhap 0 hay 1:'); disp('Thu tu cac lop mang duoc phu tren de:') Cách 1: for u=1:solop % 3 if (tam==1) VIC FAN
SWEET NOVEMBER 45 VERSION 2009 thutu0(2*u-1)=1; end thutu0(2*u)=0; x=x+1; thutu(u)=thutu0(u); end else else thutu0(2*u-1)=0; while (x chietsuat2(1)) n=1; while (x chietsuat2(1)) while (x<=solop) CAU 7 y=1; % 7 while (y<=length(buocsong)) plot(buocsong,R); csuat(2*x-1,y)=chietsuat2(y); xlabel('Buoc song (nanomet)'); csuat(2*x,y)=chietsuat1(y); ylabel('Do phan xa R'); chietsuat(x,y)=csuat(x,y); title('Do thi do phan xa R thay doi theo buoc song'); y=y+1; VIC FAN
SWEET NOVEMBER 46 VERSION 2009 Nhap do day tung lop mang (nanomet):146 KẾT QUẢ THU ĐƯỢC Nhap vao buocsongmin:1000 doday = Nhap vao buocsongmax:2000 Nhap vao buoc nhay:2 245 338 81 25 30 26 146 Danh sach cac chat co the mo phong: 1-PbTe 2-Ge 3-CdTe 4-ZnSe 5-ZnS 6-Ta2O5 Nhap vao chiet suat de:1.52 7-YbF3 8-YF3 9-SiO2 10-MgF2 11-CaF2 12- Nhap vao chiet suat moi truong:1 BaF2 Nhap vao 2 chat Chat thu nhat:6 Chat thu hai:10 Chat co chiet suat cao la chat thu: key1 = 6 Chat co chiet suat thap la chat thu: key2 = 10 Nhap vao so lop mang ban muon mo phong:7 Nhap 0 hay 1:1 Thu tu cac lop mang duoc phu tren de: solop = thutu = 7 1 0 1 0 1 0 1 Cho biet lop mang dau tien la lop cua chat co chiet suat cao hay chiet suat thap. Nhap do day tung lop mang (nanomet):245 Neu la chat co chiet suat cao thi nhap vao so 1. Nhap do day tung lop mang (nanomet):338 Neu la chat co chiet suat thap thi nhap vao so Nhap do day tung lop mang (nanomet):81 0. Nhap do day tung lop mang (nanomet):25 Nhap 0 hay 1:0 Nhap do day tung lop mang (nanomet):30 Thu tu cac lop mang duoc phu tren de: Nhap do day tung lop mang (nanomet):26 Nhap do day tung lop mang (nanomet):146 thutu = doday = 0 1 0 1 0 1 0 245 338 81 25 30 26 146 Nhap do day tung lop mang (nanomet):245 Nhap do day tung lop mang (nanomet):338 Nhap vao chiet suat de:1.52 Nhap do day tung lop mang (nanomet):81 Nhap vao chiet suat moi truong:1 Nhap do day tung lop mang (nanomet):25 Nhap do day tung lop mang (nanomet):30 Nhap do day tung lop mang (nanomet):26 VIC FAN
SWEET NOVEMBER 47 VERSION 2009 lamda=lamda+buocnhay; end buocsong % Cau 2 % Chon ra 2 chat bat ky trong so 12 chat syms PbTe Ge CdTe ZnSe ZnS Ta2O5 YbF3 YF3 SiO2 MgF2 CaF2 BaF2 chat=[PbTe Ge CdTe ZnSe ZnS Ta2O5 YbF3 YF3 SiO2 MgF2 CaF2 BaF2]; disp('Ban hay chon ra 2 chat khac nhau trong so 12 chat sau day: PbTe, Ge, CdTe, ZnSe, ZnS, Ta2O5, YbF3, YF3, SiO2, MgF2, CaF2, BaF2') MỘT SỐ BÀI MẪU THAM KHẢO disp('Chon bang cach go vao chinh xac ten % BAI LAP TRINH 3 cua chat.') % SV: Doan Quoc Huy chat1=input('Chat thu nhat ban chon la:'); % MSSV: 0513078 for a=1:length(chat) clc if (chat(a)==chat1) clear all ms1=a; % Cau 1 end % Chon vung buoc song khao sat end lamdamin=input('Nhap vao buoc song nho ms1 nhat (nm):'); chat2=input('Chat thu hai ban chon la:'); lamdamax=input('Nhap vao buoc song lon for b=1:length(chat) nhat (nm):'); if (chat(b)==chat2) while (lamdamax =lamdamin) a=[5.0781 4.0517 2.0681 2.0845 1.7951 buocnhay=input('Nhap lai buoc nhay 1.9079 0.90287 1.0527 0.76991 1.371 phai nho hon lamdamin:'); 1.1199 0.86148]; end b=[8.0749*10^7 6.0528*10^7 9.2789*10^7 disp('Mang buoc song khao sat la:') 3.4469*10^8 7.0775*10^8 6.4378*10^8 dem=1; 2.7005*10^9 4.8652*10^8 4.3456*10^9 lamda=lamdamin; 1.9465*10^8 1.1747*10^9 1.50799*10^9]; while (lamda<=lamdamax) c=[4 1 1 1 1 1 1 1 1 1 1 1]; buocsong(dem)=lamda; d=[6715 8131.5 8157.6 6715 35622 27923 xicma(dem)=1/lamda; 15717 14675 40506 132800 41379 30788]; dem=dem+1; VIC FAN
SWEET NOVEMBER 48 VERSION 2009 f=[6.0882*10^7 7.9566*10^6 1.2387*10^8 thutu1 2.7301*10^8 8.3564*10^8 2.2668*10^9 elseif (huy==thap) 1.9486*10^9 1.1942*10^9 4.4229*10^9 disp('Thu tu mang da lop cua ban la:') 2.9588*10^9 4.8958*10^9 2.3347*10^9]; dem2=1; % Chiet suat cua chat thu nhat for c2=1:solop chietsuat1=a(ms1)+b(ms1)/(c(ms1)*(xicma thutu(2*c2-1)=thap; (1)-d(ms1))^2+f(ms1)); thutu(2*c2)=cao; % Chiet suat cua chat thu hai thutu2(dem2)=thutu(dem2); chietsuat2=a(ms2)+b(ms2)/(c(ms2)*(xicma dem2=dem2+1; (1)-d(ms2))^2+f(ms2)); end % So sanh chiet suat 2 chat thutu2 if (chietsuat1>chietsuat2) else disp('Chat co chiet suat cao la:') error('Ban da go sai! Vui long go lai') chat1 end disp('Chat co chiet suat thap la:') % Cau 5 chat2 % Nhap do day cac lop else disp('Nhap vao do day (nm)') disp('Chat co chiet suat cao la:') for f1=1:solop chat2 doday(f1)=input('Do day tung lop la:'); disp('Chat co chiet suat thap la:') end chat1 doday end % Cau 6 % Cau 4 % Tao mang chiet suat thay doi theo buoc solop=input('Nhap vao so lop mang ban can song cua tung chat mo phong:') for bien=1:length(xicma) syms cao thap chon=[cao thap]; chietsuat1(bien)=a(ms1)+b(ms1)/(c(ms1)*( disp('Ban muon chon lop dau tien la chat co xicma(bien)-d(ms1))^2+f(ms1)); chiet suat cao hay thap?') disp('Neu muon la chat co chiet suat cao, chietsuat2(bien)=a(ms2)+b(ms2)/(c(ms2)*( hay go chu: cao') xicma(bien)-d(ms2))^2+f(ms2)); disp('Neu muon la chat co chiet suat thap, end hay go chu: thap') % Tao ma tran chiet suat gom co solop huy=input('Ban vui long go chu:'); dong va length(buocsong) cot if (huy==cao) if (huy==cao) disp('Thu tu mang da lop cua ban la:') if (chietsuat1(1)>chietsuat2(1)) dem1=1; for h1=1:solop for c1=1:solop n1(2*h1-1,:)=chietsuat1; thutu(2*c1-1)=cao; n1(2*h1,:)=chietsuat2; thutu(2*c1)=thap; mt_chietsuat(h1,:)=n1(h1,:); thutu1(dem1)=thutu(dem1); end dem1=dem1+1; else end for h1=1:solop VIC FAN
SWEET NOVEMBER 49 VERSION 2009 n1(2*h1-1,:)=chietsuat2; C(h3)=matrix(2,1); n1(2*h1,:)=chietsuat1; R(h3)=(abs((B(h3)- mt_chietsuat(h1,:)=n1(h1,:); C(h3))/(B(h3)+C(h3))))^2; end end end R=R*100 else % Cau 7 if (chietsuat1(1)>chietsuat2(1)) % Do thi bieu dien su phu thuoc cua do for h1=1:solop phan xa theo buoc song n1(2*h1-1,:)=chietsuat2; plot(buocsong,R); n1(2*h1,:)=chietsuat1; xlabel('Buoc song (nanomet)'); mt_chietsuat(h1,:)=n1(h1,:); ylabel('Do phan xa R (%)'); end title('Do thi bieu dien do phan xa R thay doi else theo buoc song'); for h1=1:solop grid on n1(2*h1-1,:)=chietsuat1; n1(2*h1,:)=chietsuat2; PHẦN KẾT QUẢ mt_chietsuat(h1,:)=n1(h1,:); Nhap vao buoc song nho nhat (nm):1000 end Nhap vao buoc song lon nhat (nm):2000 end Nhap vao buoc nhay (nm):2 end Ban hay chon ra 2 chat khac nhau trong so % Chiet suat de 12 chat sau day: PbTe, Ge, CdTe, ZnSe, ns=1.52; ZnS, Ta2O5, YbF3, YF3, SiO2, MgF2, % Ap dung cong thuc tinh do phan xa cua CaF2, BaF2 mang da lop thay doi theo buoc song Chon bang cach go vao chinh xac ten cua for h3=1:length(buocsong) chat. mtdv=[1 0;0 1]; Chat thu nhat ban chon la:Ta2O5 for h4=1:solop ms1 = 6 Chat thu hai ban chon la:MgF2 theta(h4)=2*pi*mt_chietsuat(h4,h3)*doday ms2 = 10 (h4)*1/(buocsong(h3)); Chat co chiet suat cao la: Q(h4)=cos(theta(h4)); chat1 = Ta2O5 Chat co chiet suat thap la: U(h4)=0+sin(theta(h4))/mt_chietsuat(h4,h3 chat2 = MgF2 )*i; Nhap vao so lop mang ban can mo phong:7 solop = 7 O(h4)=0+mt_chietsuat(h4,h3)*sin(theta(h4) Ban muon chon lop dau tien la chat co chiet )*i; suat cao hay thap? C(h4)=cos(theta(h4)); Neu muon la chat co chiet suat cao, hay go H=[Q(h4) U(h4);O(h4) C(h4)]; chu: cao mtdv=mtdv*H; Neu muon la chat co chiet suat thap, hay go end chu: thap matrix=mtdv*[1;ns]; B(h3)=matrix(1,1); VIC FAN
SWEET NOVEMBER 50 VERSION 2009 Trường hợp 1: Lớp đầu tiên là lớp có chiết suất cao Ban vui long go chu:cao Thu tu mang da lop cua ban la: thutu1 = [cao, thap, cao, thap, cao, thap, cao] Nhap vao do day (nm) Do day tung lop la:245 Do day tung lop la:338 Do day tung lop la:81 Do day tung lop la:25 Do day tung lop la:30 Do day tung lop la:26 %%% %% % BAI 3 % %% %%% Do day tung lop la:146 % NGUYEN TRUNG DUONG % MSSV: 0513061 clc clear all clear figure %%% %% % Cau 1 % %% %%% % Chon lamda_min, lamda_max, buoc_nhay lamda_min=input('Nhap buoc song min (don vi nm):'); lamda_max=input('Nhap buoc song max (don vi nm):'); buoc_nhay=input('Nhap buoc nhay (don vi nm):'); a=1; for b=lamda_min:buoc_nhay:lamda_max Trường hợp 2: Lớp đầu tiên là lớp có chiết lamda(a)=b; suất thấp xicma(a)=1/b; Ban vui long go chu:thap a=a+1; Thu tu mang da lop cua ban la: end thutu1 = lamda [thap, cao, thap, cao, thap, cao, thap] %%% %% % Cau 2 % %% %%% Nhap vao do day (nm) % Chon chat: chon ra 2 chat bat ky trong so 12 Do day tung lop la:245 chat Do day tung lop la:338 disp('Ma so cac chat mo phong:') disp('So 1: PbTe, So 2: Ge, So 3: CdTe, So 4: Do day tung lop la:81 ZnSe') Do day tung lop la:25 disp('So 5: ZnS, So 6: Ta2O5, So 7: YbF3, So Do day tung lop la:30 8: YF3') Do day tung lop la:26 disp('So 9: SiO2, So 10: MgF2, So 11: CaF2, Do day tung lop la:146 So 12: BaF2') chat_1=input('Nhap vao chat thu nhat:'); while (chat_1 12) VIC FAN
SWEET NOVEMBER 51 VERSION 2009 chat_1=input('Sai ma so!. Vui long nhap lai') chat_1 end tam=n1; chat_2=input('Nhap vao chat thu hai:'); n1=n2; while (chat_2 12) n2=tam; chat_2=input('Sai ma so!. Vui long nhap lai') end end %%% %% % Cau 4 % %% %%% %%% %% % Cau 3 % %% %%% disp('Neu chon lop dau tien la lop co chiet suat % Mang du lieu de tinh chiet suat cua 12 chat cao, nhap so 1') a1=[5.0781 4.0517 2.0681 2.0845 1.7951 disp('Neu chon lop dau tien la lop co chiet suat 1.9079 0.90287 1.0527 0.76991 1.371 1.1199 thap, nhap so 0') 0.86148]; nhap=input('Nhap:'); a2=[8.0749*10^7 6.0528*10^7 9.2789*10^7 while (nhap~=0)&(nhap~=1) 3.4469*10^8 7.0775*10^8 6.4378*10^8 nhap=input('Chi nhap 0 hoac 1. Vui long 2.7005*10^9 4.8652*10^8 4.3456*10^9 nhap lai:'); 1.9465*10^8 1.1747*10^9 1.50799*10^9]; end a3=[4 1 1 1 1 1 1 1 1 1 1 1]; so_lop=input('Nhap vao so lop mang:'); a4=[6715 8131.5 8157.6 6715 35622 27923 disp('Thu tu cac lop mang la:') 15717 14675 40506 132800 41379 30788]; if (nhap==0) a5=[6.0882*10^7 7.9566*10^6 1.2387*10^8 for d1=1:2:so_lop 2.7301*10^8 8.3564*10^8 2.2668*10^9 mang(d1)=0; 1.9486*10^9 1.1942*10^9 4.4229*10^9 mang(d1+1)=1; 2.9588*10^9 4.8958*10^9 2.3347*10^9]; end for a6=1:length(xicma) else mau1(a6)=a3(chat_1)*(xicma(a6)- for d2=1:2:so_lop a4(chat_1))^2+a5(chat_1); mang(d2)=1; n1(a6)=a1(chat_1)+a2(chat_1)/mau1(a6); mang(d2+1)=0; mau2(a6)=a3(chat_2)*(xicma(a6)- end a4(chat_2))^2+a5(chat_2); end n2(a6)=a1(chat_2)+a2(chat_2)/mau2(a6); for d3=1:so_lop end thutu(d3)=mang(d3); disp('Chiet suat cua chat thu nhat la:') end n1 thutu disp('Chiet suat cua chat thu hai la:') %%% %% % Cau 5 % %% %%% n2 disp('Do day cac lop lan luot la:') disp('So sanh 2 chiet suat, ta ket luan:') for bien=1:so_lop if (n1(1)>n2(1)) do_day(bien)=input('Nhap:'); disp('Chat co chiet suat cao la chat co ma end so:') %%% %% % Cau 6 % %% %%% chat_1 % Ket hop 2 mang chiet suat n1 va n2 thanh disp('Chat co chiet suat thap la chat co ma mang n duy nhat so:') if (nhap==1) chat_2 for t1=1:2:so_lop else for t2=1:length(lamda) disp('Chat co chiet suat cao la:') c_s(t1,t2)=n1(t2); chat_2 c_s(t1+1,t2)=n2(t2); disp('Chat co chiet suat thap la:') end VIC FAN
SWEET NOVEMBER 52 VERSION 2009 end Nhap buoc nhay (don vi nm):2 else Ma so cac chat mo phong: for t1=1:2:so_lop So 1: PbTe, So 2: Ge, So 3: CdTe, So 4: ZnSe for t2=1:length(lamda) So 5: ZnS, So 6: Ta2O5, So 7: YbF3, So 8: c_s(t1,t2)=n2(t2); YF3 c_s(t1+1,t2)=n1(t2); So 9: SiO2, So 10: MgF2, So 11: CaF2, So 12: end BaF2 end Nhap vao chat thu nhat:6 end Nhap vao chat thu hai:10 for t3=1:so_lop So sanh 2 chiet suat, ta ket luan: for t4=1:length(lamda) Chat co chiet suat cao la chat co ma so: n(t3,t4)=c_s(t3,t4); chat_1 = end 6 end Chat co chiet suat thap la chat co ma so: csd=1.52 % Chiet suat de chat_2 = for tt=1:length(lamda) 10 A=[1 0;0 1]; Neu chon lop dau tien la lop co chiet suat cao, for t=1:so_lop nhap so 1 Neu chon lop dau tien la lop co chiet suat thap, theta(t)=2*pi*n(t,tt)*do_day(t)*1/(lamda(tt)); nhap so 0 A1(t)=cos(theta(t)); A2(t)=0+sin(theta(t))/n(t,tt)*i; Nhap:1 A3(t)=0+n(t,tt)*sin(theta(t))*i; Nhap vao so lop mang:7 A4(t)=cos(theta(t)); Thu tu cac lop mang la: A0=[A1(t) A2(t);A3(t) A4(t)]; thutu = A=A*A0; 1 0 1 0 1 0 1 end Do day cac lop lan luot la: A5=A*[1;csd]; Nhap:245 B(tt)=A5(1,1); Nhap:338 C(tt)=A5(2,1); Nhap:81 A6(tt)=B(tt)-C(tt); Nhap:25 A7(tt)=B(tt)+C(tt); Nhap:30 R(tt)=(abs(A6(tt)/A7(tt)))^2; Nhap:26 end Nhap:146 R %%% %% % Cau 7 % %% %%% plot(lamda,R); xlabel('Lamda (nm)'); ylabel('Do phan xa R'); title('Do thi do phan xa R thay doi theo buoc song'); grid on PHẦN BÁO CÁO Nhap buoc song min (don vi nm):1000 Nhap buoc song max (don vi nm):2000 VIC FAN
SWEET NOVEMBER 53 VERSION 2009 disp('DANH SACH 12 CHAT CO THE MO Nhap:0 PHONG VA MA SO TUONG UNG:') Nhap vao so lop mang:7 disp('PbTe Ge CdTe ZnSe ZnS Ta2O5 YbF3 Thu tu cac lop mang la: YF3 SiO2 MgF2 CaF2 BaF2') thutu = disp(' 1 2 3 4 5 6 7 8 9 10 11 0 1 0 1 0 1 0 12') Do day cac lop lan luot la: disp('MOI BAN CHON RA 2 CHAT BAT KY Nhap:245 BANG CACH NHAP VAO MA SO TUONG Nhap:338 UNG') Nhap:81 disp('CHAT THU NHAT LA:') Nhap:25 NUM1=input(''); Nhap:30 disp('CHAT THU HAI LA:') Nhap:26 NUM2=input(''); Nhap:146 % CAU 3 K1=[5.0781 4.0517 2.0681 2.0845 1.7951 1.9079 0.90287 1.0527 0.76991 1.371 1.1199 0.86148]; K2=[8.0749*10^7 6.0528*10^7 9.2789*10^7 3.4469*10^8 7.0775*10^8 6.4378*10^8 2.7005*10^9 4.8652*10^8 4.3456*10^9 1.9465*10^8 1.1747*10^9 1.50799*10^9]; K3=[4 1 1 1 1 1 1 1 1 1 1 1]; K4=[6715 8131.5 8157.6 6715 35622 27923 15717 14675 40506 132800 41379 30788]; K5=[6.0882*10^7 7.9566*10^6 1.2387*10^8 2.7301*10^8 8.3564*10^8 2.2668*10^9 % HO TEN SV: PHAM DANG KHOA 1.9486*10^9 1.1942*10^9 4.4229*10^9 MSSV:0513009 LOP: 2.9588*10^9 4.8958*10^9 2.3347*10^9]; 05VLUD2 for RUN=1:length(BS) % MOM1(RUN)=K3(NUM1)*(1/BS(RUN)- K4(NUM1))^2+K5(NUM1); clc clear all INDEX1(RUN)=K1(NUM1)+K2(NUM1)/MO % CAU 1 M1(RUN); end BSmin=input('NHAP BUOC SONG NHO IND1=mean(INDEX1) NHAT (NANOMET):'); for RUN=1:length(BS) BSmax=input('NHAP BUOC SONG LON MOM2(RUN)=K3(NUM2)*(1/BS(RUN)- NHAT (NANOMET):'); K4(NUM2))^2+K5(NUM2); STEP=input('NHAP BUOC NHAY (NANOMET):'); INDEX2(RUN)=K1(NUM2)+K2(NUM2)/MO BS=[BSmin:STEP:BSmax] M2(RUN); % CAU 2 end IND2=mean(INDEX2) if (IND1>IND2) VIC FAN
SWEET NOVEMBER 54 VERSION 2009 disp('CHAT CO CHIET SUAT CAO LA disp('LOW') CHAT:');NUM1 end disp('CHAT CO CHIET SUAT THAP LA disp('HIGH') CHAT:');NUM2 end else else disp('CHAT CO CHIET SUAT CAO LA error('TEST CHI BANG 1 HOAC 2. VUI CHAT:');NUM2 LONG CHAY LAI TU DAU.') disp('CHAT CO CHIET SUAT THAP LA end CHAT:');NUM1 disp(' ') CHANGE=INDEX1 % CAU 5 INDEX1=INDEX2 INDEX2=CHANGE disp('HAY NHAP VAO DO DAY MOI LOP') end ONE=1; % CAU 4 while (ONE<=LAYER) DODAY(ONE)=input('DO DAY LA:'); disp('BAN CHON LOP DAU TIEN LA LOP ONE=ONE+1; CO CHIET SUAT THAP HAY CAO?') end disp('NEU THAP THI NHAP TEST = 1') DODAY disp('NEU CAO THI NHAP TEST = 2') % CAU 6 TEST=input('TEST = '); LAYER=input('NHAP VAO SO LOP MANG S=1.52; BAN MUON MO PHONG:'); if (TEST==1) disp('THU TU PHU MANG LA:') for S1=1:length(BS) disp(' ') MASK=eye(2,2); if (TEST==1) for S2=1:LAYER if (mod(LAYER,2)==0) if (mod(S2,2)~=0) for MOVE=1:LAYER/2 disp('LOW') THETA(S2)=2*pi*INDEX2(S1)*DODAY(S2) disp('HIGH') *1/(BS(S1)); end TITAN1(S2)=cos(THETA(S2)); else for MOVE=1:LAYER/2 TITAN2(S2)=0+sin(THETA(S2))/INDEX2(S1 disp('LOW') )*i; disp('HIGH') end TITAN3(S2)=0+INDEX2(S1)*sin(THETA(S2) disp('LOW') )*i; end TITAN4(S2)=cos(THETA(S2)); elseif (TEST==2) TITAN=[TITAN1(S2) if (mod(LAYER,2)==0) TITAN2(S2);TITAN3(S2) TITAN4(S2)]; for MOVE=1:LAYER/2 MASK=MASK*TITAN; disp('HIGH') else disp('LOW') end THETA(S2)=2*pi*INDEX1(S1)*DODAY(S2) else *1/(BS(S1)); for MOVE=1:LAYER/2 TITAN1(S2)=cos(THETA(S2)); disp('HIGH') VIC FAN
SWEET NOVEMBER 55 VERSION 2009 TITAN2(S2)=0+sin(THETA(S2))/INDEX1(S1 TITAN3(S2)=0+INDEX2(S1)*sin(THETA(S2) )*i; )*i; TITAN4(S2)=cos(THETA(S2)); TITAN3(S2)=0+INDEX1(S1)*sin(THETA(S2) TITAN=[TITAN1(S2) )*i; TITAN2(S2);TITAN3(S2) TITAN4(S2)]; TITAN4(S2)=cos(THETA(S2)); MASK=MASK*TITAN; TITAN=[TITAN1(S2) end TITAN2(S2);TITAN3(S2) TITAN4(S2)]; end MASK=MASK*TITAN; MODERN=MASK*[1;S]; end B(S1)=MODERN(1,1); end C(S1)=MODERN(2,1); MODERN=MASK*[1;S]; R(S1)=(abs((B(S1)- B(S1)=MODERN(1,1); C(S1))/(B(S1)+C(S1))))^2; C(S1)=MODERN(2,1); end R(S1)=(abs((B(S1)- end C(S1))/(B(S1)+C(S1))))^2; % CAU 7 end else plot(BS,R); for S1=1:length(BS) xlabel('BUOC SONG (NANOMET)'); MASK=eye(2,2); ylabel('DO PHAN XA R'); for S2=1:LAYER title('SU PHU THUOC CUA R THEO BUOC if (mod(S2,2)~=0) SONG'); grid on THETA(S2)=2*pi*INDEX1(S1)*DODAY(S2) *1/(BS(S1)); TITAN1(S2)=cos(THETA(S2)); KET QUA THU DUOC TITAN2(S2)=0+sin(THETA(S2))/INDEX1(S1 NHAP BUOC SONG NHO NHAT )*i; (NANOMET):1000 NHAP BUOC SONG LON NHAT TITAN3(S2)=0+INDEX1(S1)*sin(THETA(S2) (NANOMET):2000 )*i; NHAP BUOC NHAY (NANOMET):2 TITAN4(S2)=cos(THETA(S2)); DANH SACH 12 CHAT CO THE MO TITAN=[TITAN1(S2) PHONG VA MA SO TUONG UNG: TITAN2(S2);TITAN3(S2) TITAN4(S2)]; PbTe Ge CdTe ZnSe ZnS Ta2O5 YbF3 YF3 MASK=MASK*TITAN; SiO2 MgF2 CaF2 BaF2 else 1 2 3 4 5 6 7 8 9 10 11 12 THETA(S2)=2*pi*INDEX2(S1)*DODAY(S2) MOI BAN CHON RA 2 CHAT BAT KY *1/(BS(S1)); BANG CACH NHAP VAO MA SO TUONG TITAN1(S2)=cos(THETA(S2)); UNG CHAT THU NHAT LA: TITAN2(S2)=0+sin(THETA(S2))/INDEX2(S1 6 )*i; CHAT THU HAI LA: 10 VIC FAN
SWEET NOVEMBER 56 VERSION 2009 IND1 = 2.1192 IND2 = 1.3805 CHAT CO CHIET SUAT CAO LA CHAT: NUM1 = 6 CHAT CO CHIET SUAT THAP LA CHAT: NUM2 = 10 BAN CHON LOP DAU TIEN LA LOP CO CHIET SUAT THAP HAY CAO? NEU THAP THI NHAP TEST = 1 TEST = 1 NEU CAO THI NHAP TEST = 2 NHAP VAO SO LOP MANG BAN MUON TEST = 2 MO PHONG:7 NHAP VAO SO LOP MANG BAN MUON THU TU PHU MANG LA: MO PHONG:7 THU TU PHU MANG LA: LOW HIGH HIGH LOW LOW HIGH HIGH LOW LOW HIGH HIGH LOW LOW HIGH HAY NHAP VAO DO DAY MOI LOP DO DAY LA:245 HAY NHAP VAO DO DAY MOI LOP DO DAY LA:338 DO DAY LA:245 DO DAY LA:81 DO DAY LA:338 DO DAY LA:25 DO DAY LA:81 DO DAY LA:30 DO DAY LA:25 DO DAY LA:26 DO DAY LA:30 DO DAY LA:146 DO DAY LA:26 DODAY = DO DAY LA:146 245 338 81 25 30 26 146 DODAY = 245 338 81 25 30 26 146 VIC FAN
SWEET NOVEMBER 57 VERSION 2009 chiết suất đế là 1.52. Chiết suất môi trường là 1. 3) Thực hiện thuật giải N-Square Scan để tìm độ dày từng lớp màng sao cho độ phản xạ đạt được gần bằng giá trị Rtarget nhất. Độ phản xạ R của màng theo bước sóng dựa vào công thức: B C 2 R  B C Trong đó: q isin j B cos j 1 n j C  n BÀI LẬP TRÌNH 5 (Cảnh báo! Đây là bài j 1 in sin cos s KHÓ nhất trong loạt bài mô phỏng quang học j j j bằng Matlab! Lớp 04 tác giả ngày xưa không 2 n jd j Với:  ai làm nổi!!! Ngoài ra, thời gian để chạy bài m  này trung bình từ 30 phút đến 2 tiếng, một cách Hàm so sánh giá trị: tốt để kiểm tra sức mạnh máy tính ) n 1/ 2 1 2  Các yêu cầu trong bài lập trình: F k  Rt arg et R k   1) Thực hiện lại từ yêu cầu 1 đến yêu cầu 4 của n k 1  bài lập trình 4. Các yêu cầu cụ thể như sau: n là số phần tử trong mảng bước sóng a) Lập trình chọn vùng bước sóng khảo sát Hướng dẫn về thuật giải như sau: (Ví dụ đơn (lamdamin, lamdamax, bước nhảy) giản với 3 lớp màng) b) Lập trình chọn ra 2 chất bất kỳ trong số Lượt chạy 1: 12 chất: PbTe, Ge, CdTe, ZnSe, ZnS, Ta2O5, - Cho lớp 1 chạy (từ 2 → 300nm), cố định lớp YbF3, YF3, SiO2, MgF2, CaF2, BaF2. 2 và lớp 3 ở 2nm. Tìm được độ dày d11 của lớp In kết quả với 2 chất là Ta2O5 và MgF2 1 tại đó Fmin. c) So sánh chiết suất của 2 chất vừa chọn ở - Giữ lớp 1 bằng d11, cho lớp 2 chạy (từ 2 → câu 2 tại phần tử đầu tiên của mảng bước sóng. 300nm), cố định lớp 3 ở 2nm. Tìm được độ dày Chất nào có chiết suất lớn hơn sẽ là chất chiết d21 của lớp 2 tại đó Fmin. suất cao, chất nào có chiết suất nhỏ hơn sẽ là - Giữ lớp 1 bằng d11, lớp 2 bằng d21, cho lớp 3 chất chiết suất thấp. chạy (từ 2 → 300nm), tìm được độ dày d31 của d) Lập trình chọn chất đầu tiên là chất chiết lớp 3 tại đó Fmin suất cao hay chất chiết suất thấp. Chọn số lớp Như vậy, 3 độ dày tối ưu của lượt chạy 1 là: màng cần mô phỏng. Sắp xếp xen kẽ các lớp d11, d21, d31. (Số 1 ở cuối thể hiện lượt chạy màng theo thứ tự chiết suất cao – chiết suất 1) thấp – chiết suất cao – hay chiết suất thấp – Lượt chạy 2: chiết suất cao – chiết suất thấp – - Giữ lớp 1 bằng d11, lớp 3 bằng d31, cho lớp 2 In kết quả với số lớp là 7. chạy (từ 2 → 300nm). Tìm được độ dày d22 2) Nhập vào độ phản xạ Rtarget mà bạn mong của lớp 2 tại đó Fmin. muốn màng đa lớp của bạn đạt được. Nhập vào - Giữ lớp 2 bằng d22, giữ lớp 1 bằng d11, chạy khoảng giới hạn độ dày (độ dày nhỏ nhất, độ lớp 3 (từ 2 → 300nm). Tìm được độ dày d32 dày lớn nhất). của lớp 3 tại đó Fmin. Chạy chương trình với Rtarget = 0.005. - Giữ lớp 2 bằng d22, giữ lớp 3 bằng d32, cho Khoảng giới hạn độ dày: 2 → 300nm. Cho biết lớp 1 chạy (từ 2 → 300nm). Tìm được độ dày d12 của lớp 1 tại đó Fmin. VIC FAN
SWEET NOVEMBER 58 VERSION 2009 Như vậy, 3 độ dày tối ưu của lượt chạy 2 là: d12, d22, d32. Lượt chạy 3: - Giữ lớp 1 bằng d12, giữ lớp 2 bằng d22, cho lớp 3 chạy (từ 2 → 300nm). Tìm được độ dày d33 của lớp 3 tại đó Fmin. - Giữ lớp 2 bằng d22, giữ lớp 3 bằng d33, cho lớp 1 chạy (từ 2 → 300nm). Tìm được độ dày d13 của lớp 1 tại đó Fmin. - Giữ lớp 3 bằng d33, giữ lớp 1 bằng d13, cho lớp 2 chạy (từ 2 → 300nm). Tìm được độ dày d23 của lớp 2 tại đó Fmin. Như vậy, 3 độ dày tối ưu của lượt chạy 2 là: d13, d23, d33. 3 độ dày tối ưu cần tìm của cả bài toán: d13, d23, d33. 4) Từ 7 độ dày tối ưu của 7 lớp màng tìm được, vẽ đồ thị độ phản xạ R thay đổi theo bước sóng. Hướng dẫn: CÂU 1: Các bạn copy và paste lại từ bài lập trình 3 ở trên. (Câu này dễ, hy vọng bạn nào cũng có thể làm được) CÂU 2: Câu này cũng dễ luôn! Bạn đọc có thể tham khảo đoạn chương trình bên dưới VIC FAN