-
Notifications
You must be signed in to change notification settings - Fork 1
/
blueredyellow.m
88 lines (67 loc) · 2.11 KB
/
blueredyellow.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
function [colors] = blueredyellow(K)
% High contrast colormap that prints to grayscale correctly
% ======================================================================
% Copyright (c) 2012 David Weiss
%
% Permission is hereby granted, free of charge, to any person obtaining
% a copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to
% permit persons to whom the Software is furnished to do so, subject to
% the following conditions:
%
% The above copyright notice and this permission notice shall be
% included in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
% EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
% NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
% LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
% OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
% WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
% ======================================================================
if nargin==0
K = 64;
end
MR=[0,0;
0.1,0;
% 0.3,0.5;
0.7,1;
% 0.8,0.8;
1,1];
MG=[0,0;
0.3,0;
0.8,0.7;
1,1];
MB=[0,0.1;
0.3,0.6;
% 0.5,0.3;
0.7,0.2;
1, 0.5];
colors = colormapRGB(K, MR,MG,MB);
return;
%%
x = linspace(-20,20,500);
%compose a grid
[xx,yy] = meshgrid(x,x);
% calculate the radius r for xx and yy pairs
% as the airy disc has a circular symmetry
r = sqrt( xx.^2 + yy.^2 );
%calculate the electrical field
AiryField = besselj(1,r)./r;
% intensity = abs( field ) ^2
AiryIntensity = abs(AiryField).^2;
% plot them
figure(1)
imagesc(x,x,AiryField)
colormap(blueredyellow);
figure(2)
imagesc(x,x,AiryField)
colormap(jet);
%%
imagesc(rand(100,100)>0.5);
%%
% figure(2)
% imagesc(x,x,AiryIntensity);