• 沒有找到結果。

4. Pupil filters designed for simultaneously achieving super-resolution for

4.4 Set-up of The Four-zone Pupil Filter

4.5.2 Simulation Verification

Step 1:

With the condition a2+b2=1, we adjust the phase factors and the radii of these four zones. For CD surface, t1 and t3 are still assigned both to be 1; and φ1 and φ3 are assigned to be 0 and 0.4π. For DVD surface, T2 and T4 are also assigned both to be 1; and Ψ2 and Ψ4 are assigned to be 0 and 0.35π. The value of radius c is first assumed to be 0.7, while the value of radius a varies from 0.4 to 0.6, and the value of radius b is obtained from the relation b= (1-a2)1/2.

With all the factors being settled down, we obtain the relation among the radius a and the transverse gains and Strehl ratio of such a system for both CD and DVD cases, as shown in Fig. 4.7(a) and 4.7(b), respectively. It can be observed that with increase in radius a, the transverse gain decreases for CD, while one increases first and then decreases for DVD. In the range that a∈[0.4, 0.6], the gains are all greater than 1, which means that superresolution property can be realized simultaneously for both CD and DVD. It can also be seen that, for the case of CD, the Strehl ratio increases with the increment of the radius a, but decreases for the case of DVD.

Fig. 4.7: 1) Relation among the transverse gain and the radius of the first zone.

2) Relation among the Strehl ratio and the radius of the first zone.

(a) for the case of CD (b) for the case of DVD

Figures 4.8(a) and 4.8(b), respectively show the transverse intensity distributions for three particular solutions for CD and DVD, corresponding to three different kinds of set-up of the pupil filter, in comparison to the case of clear pupil. The intensity has been normalized to the clear pupil size.

It should be carefully minded that the simulated results of the intensity distribution, shown in Fig. 4.8, disagree with those derived from Eqs. 4.11 and 4.12.

For instance, when assuming the value of radius a to be 0.6, the predicted value of the Strehl ratio for CD derived from Eq. 4.11 is about 0.9, while that for DVD is approximately equal to 0.5.

Fig. 4.8: The transverse intensity distributions 1) with the clear pupil, 2) with the designed filter, for a=0.5 (solid curve), a=0.57 (dashed curve), a=0.6 (dotted curve) (a) for the case of CD (b) for the case of DVD.

However, as shown in figure 4.8, the central peak value of the normalized intensity distribution for CD is, in fact, less than 0.4, while that for DVD is less than

0.5. Besides, for CD, the set-up of the filter doesn’t achieve the goal of superresolution, but even enlarges the spot size in the focal plane. For DVD, though transverse superresolution is obtained, the transverse side-lobe is tremendously worsened in the focal plane. It is noticed that, for both cases of CD and DVD, the contrast is worsened, which makes it much more difficult in the practical application of reading the data from the disk.

(a) (b)

Fig. 4.9: Relation among the displacement of focus in the axial direction and the radius of the first zone. (a) for the case of CD (b) for the case of DVD

It can be clearly be seen in Fig. 4.9 that the computed value of the axial displacement, uF, is in fact quite apart from zero for both cases of CD and DVD.

The place where maxima intensity occurs has been shifted away from the geometric focus, so that the derived forms of Strehl ratio and transverse gain, based on 2nd-order approximation, become incorrect in describing the focal behavior for the designed cases here.

Step 2:

Now, the intensity distributions along the transverse direction for CD and DVD

at the paraxial focus are obtained directly from the computation of the diffraction theory, shown in Fig. 4.10(a) and (b), respectively. We can see that the results correspond to what we have shown in Fig. 4.6. For CD, the goal of superresolution isn’t achieved, but even enlarges the spot size; for DVD, though transverse superresolution is obtained, the transverse side-lobe is worsened.

Fig. 4.10 (a)

Fig. 4.10(b)

Fig. 4.10: The transverse intensity distributions 1) with the clear pupil, 2) with the designed filter, for a=0.6 (a) for the case of CD (b) for the case of DVD.

We are now interested in where the best image plane (BIP) appears after we add the designed filter into the system. Radial intensity scans at various axial coordinates are shown in Fig. 4.11(a) and(b), for CD and DVD, respectively. The range of axial coordinates, u, explored here is -5 to 5, which would certainly seem to cover the transition region we are interested in. From the figure, it seems reasonable to state that the BIP, where best image performance is gotten for both CD and DVD, occurs around u= 3.

For CD, we can see that the peak value of the main lobe is 0.5 and the side-lobe is extremely small when the radial distance, v, ranges between -10 to 10. The size of blur circle remains the same as the one without adding the filter. For DVD,

transverse superresolution is obtained, but the increased transverse side-lobe will still be a big concern. The contrast of the system performance will be lowered down.

Fig. 4.11: Radial intensity scans at various planes (the range of axial coordinates, u, explored here is -6 to 6) (a) for the case of CD (b) for the case of DVD

Step 3:

When dealing with the case in which the best image plane is not near the paraxial focus, it’ll be more proper to use the modified 2nd order approximation method. For CD, we find that the BIP has been shifted to u=3.527, while that occurs at u=3.1 for DVD. Once the location of the BIP is found, we then fine tune the transmittances of the zones. Here, in this case, we set the transmittance of first zone for CD surface to be 0.7. The computed results of the superresolution factors are shown in Fig. 4.12.

(a) (b)

Fig. 4.12: 1) Relation among the displacement and the radius of the first zone.

2) Relation among the transverse gain and the radius of the first zone.

3) Relation among the Strehl ratio and the radius of the first zone.

(a) for the case of CD (b) for the case of DVD

It can be seen that the computed value of the axial displacement is very close to zero when the radius of the first zone is set to be 0.6. The values of transverse gains are all greater than 1 for CD and DVD, which means that superresolution property can be achieved in this case. However, as a drawback, the Strehl ratios for both CD and DVD have dropped greatly, lower than half value of that without a filter.

Step 4:

Similar to what we have done in step 2, now we try to verify the accuracy of the computed results gotten in step 3. The intensity distributions along the transverse direction for CD and DVD at the shifted focus are shown in Fig. 4.13(a) and (b), respectively.

Fig. 4.13 The transverse intensity distributions 1) with the clear pupil, 2) with the designed filter, for radius a=0.5 ; transmittance t1= 0.7 (a) for the case of CD (b) for the case of DVD.

For CD, we find that the transverse gain, in fact, is equal to 1.038, being lower than the expected value; while for DVD, that is equal to 1.259. Besides, the size of the main-lobe of the diffracted pattern has been narrowed down, achieving superresolution in both cases.

Within the range , we can see that, for CD, the transverse side-lobes become relatively smaller than the central peak value (being approximately 1% of the central peak value). But, for DVD, the transverse side-lobe isn’t alleviated but enhanced, which leads to worsen the contrast of the final image.

] 10 , 10 [ '∈ −

u

4.6 Conclusions

In this section, we have shown that, with the use of a rotationally symmetric four-zone pupil filter, transverse superresolution can be both realized for two different wavelengths. Notice that the expressions of the factors, like gains, Strehl ratio, and axial displacement, have to be modified for the case in which the best image plane is not near the paraxial focus. It should also be mentioned that, the transverse sidelobe is somehow troublesome, especially for the case of DVD. In the future work, the main aim in optical superresolution is to reduce the main-lobe size of the point-spread function while increasing the central intensity and suppressing the sidelobes.

Chapter 5

5.1 Conclusions

In the first part of this thesis, we have first presented a radially symmetric phase-only filter to help alleviate the effects caused by the fluctuation of third- and fifth-order spherical aberrations simultaneously. It can be clearly seen that the system tolerance to SA5 is improved by several times. Meanwhile, the trade-off issue, i.e., a reduction in intensity and the MTF, has also been discussed in details. These reductions form the baselines in practical applications of the use of phase filter.

Following the ideas mentioned above, we get interested in seeking a way to enhance the tolerance of the system to defocus or SA3 for both wavelengths. Phase pupil filters used to minimize the variation of Strehl ratio with defocus and third-order spherical aberration (SA3) for dual wavelengths are being discussed. Notice that if the amount of variation of the frequency is small, we may find a form of pupil phase function to achieve that goal.

In the last part of this thesis, we have presented a four-zone filter to help achieve transverse superresolution for two different wavelengths, simultaneously. There still

remains some room for improvement. While reducing the main-lobe size of the point-spread function, we need to search for a way to increase the central intensity and suppress the sidelobes.

5.2 Future Work

For the work in Chapter 2:

It’s interesting that if we can derive the form of the pupil phase function that extended to even higher order spherical aberration (say, the seventh order spherical aberration, SA7) by using the method we mentioned here. Besides, an investigation similar to that above may also be carried out for the odd-aberration case, i.e., the focal shift, primary and secondary circular coma.

For the work in Chapter 3:

It may be worthwhile of considering the use of hybrid surface lens, an objective lens of combining aspheric surfaces of DVD and CD. In that way, for those two different frequencies, the tolerance of the axial intensity to defocus or SA3 may be simultaneously enhanced.

For the work in Chapter 4:

A global optimization process is needed in the design of the superresolving pupil plate. Once we develop up a reliable optimization method, it’ll be much easier for us to find out the best set-up configuration to meet the requirement.

Appendix A: Derivation of The Pupil Phase Functions

A.1 Derivation of the phase filter that has been developed for W040=W020=0 By substituting of Eq.(2.9) into Eq.(2.6), it leads to the expression:

0

Note that Eq. (A 1.1) is equivalent to the following equation:

0 The solution of the above equation will be:

]}

By choosing α=-5, C = -3/8, and C1=0, a simple solution form will be obtained, which is expressed as:

Note that Eq.(A 1.7) is a particular solution of Eq.(A 1.6), being one of the possible pupil phase functions. With a change of variables, we have the phase pupil function as:

respectively in Eq.(B 1.7)

A.2 Derivation of the phase filter with Maréchal treatment By substituting Eq.(2.14) into Eq.(2.6) leads to the expression:

) ] 6 } 0

Note that Eq. (A 2.1) is equivalent to the following equation:

0 The solution of the above equation will be:

) ] }

Notice that . With a change of variables, the pupil phase function is expressed as the following form:

2

Appendix B: More Information for 2nd order Approximation Theory B.1 Further Discussion to The 2nd-order Approximation Theory

We now reexamine the derivation of Eqs. 4.10 to 4.13. As a reminder, in the focal plane, the field distribution can be expressed as [7]:

=

1 (B 1.1)

Notice that the Bessel function of the first kind is expressed as:

For small distances from the focus, we can expand the expressions for the focal-plane and axial amplitudes as a power series:

Omitting the higher order terms, we get:

0 2 1 4 2

So the transverse (focal) plane intensity will be:

0 2 Re( 0 1 ) 2 Similarly, the axial intensity will be:

0 2 0 1 [Re( 2 0 | 1|2)] 2

B.2 Modification to The Expressions of The Axial and The Transverse Gain

In the recent study, Silvia Ledesma and Juan Campos have extended the expressions for the axial and the transverse gain to the case in which the best image plane is not near the paraxial focus [20]. The content of the study is shown below.

For the case in which the best image plane is not near the paraxial focus, the expressions for the axial gain, transverse gain, and the Strehl ratio need some modifications. Recall that the field along the axis is:

=

1 (B 2.1)

Where u is the axial coordinate centered at the focal plane without the filter. By evaluating |U(0,u)2| numerically from Eq. B 2.1, we find the position umax where the axial intensity is maximum. Then we calculate those factors of superresolution by use of the expansions around this point.

The second-order expansion for the axial response around umax will be:

1 + (B 2.2)

The nth moments of the pupil around umax is defined as:

=

1 (B 2.3)

Now the terms taken into account is just up to second order in u’= u- umax, Then the axial intensity is approximated as:

[| |' Re( ' '*) ' ]

For transverse response, we expand the field to second order corresponding to umax:

=

1

Then the transverse intensity can be expressed as

Re( '* ') ]

Note that u0 corresponds to the center of the parabola defined in Eq. 4.41:

* 1 2

Since u0 is measured from the BIP centered at umax, so its values will be close to zero for most functions of an optical system. Thus, the superresulution factors around umax result in:

Appendix C: MTLAB Source Codes

C.1: The Strehl ratio as a function of aberration coefficient W060

% --- TOPIC: Strehl ratio, S(w20,w40,w60), v.s. W060 --- %

% --- Drafted by Chih-Yun Chan, 2nd version, test ok --- %

% --- Used in Section 2 --- %

u1 = zeros(81,1); % Field u1 -- with ideal aperture u2 = zeros(81,1); % Field u2 – with the designed filter I1 = zeros(81,1); % Intensity I1 -- with ideal aperture I2 = zeros(81,1); % Intensity I2 -- with ideal aperture B = 7*pi; % assign parameters

Bo = 0.3*7*pi;

% implementing the integral for i = 0 : 1 : 80

for x = 0 : 0.001 : 1

y1 = x*exp(j*2*pi*(((i./10)-4)*(x^6)));

y2 = x*exp(j*2*pi*(((i./10)-4)*(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u1(i+1)= u1(i+1)+ 0.001*y1;

u2(i+1)= u2(i+1)+ 0.001*y2;

end end

I1 = u1.*conj(u1)*4*pi*pi/10;

I1 = I1/max(I1); % Normalize the intensity I2 = u2.*conj(u2)*4*pi*pi/10;

I2 = I2/max(I2); % Normalize the intensity

k = -4 : 0.1 : 4;

plot(k,I1,k,I2)

Appendix C: MATLAB Source Codes (CONTINUED) C.2: Plot of the shape of the designed pupil phase function

% --- TOPIC: The shape of the designed pupil phase function --- %

% --- Drafted by Chih-Yun Chan, 1st version, test ok --- %

% --- Used in Section 2 --- %

Th = zeros(1001,1); %Thickness as a function of the radius B = 5*pi;

Bo = 0.3*B;

for m = 1:1:1001

Th(m,1)=(Bo*((((m-501)/500)^2)^3)+B*((((m-501)/500)^2)^3)*log(((m-501)/500) +1e-8));

end

k=-1:2/1000:1 % Normalized radial coordinate

plot(k,Th);

Appendix C: MATLAB Source Codes (CONTINUED) C.3: Computed MTF with initial setting W020=0 and W040=0

% --- TOPIC: Computed MTF with initial setting W020=0 and W040=0 --- %

% --- Drafted by Chih-Yun Chan, 2nd version, test ok --- %

% --- Used in Section 2 --- %

Q1= zeros(201,201); % Field Q1 -- with w060=0 lambda Q2= zeros(201,201); % Field Q2 -- with w060=0.5 lambda Q3= zeros(201,201); % Field Q3 -- with w060=1 lambda Q4= zeros(201,201); % Field Q4 -- with w060=0 lambda mtf1 = zeros(80,1); % MTF1 -- with w060=0 lambda mtf2 = zeros(80,1); % MTF1 -- with w060=0 lambda mtf3 = zeros(80,1); % MTF1 -- with w060=0 lambda mtf4 = zeros(80,1); % MTF1 -- with w060=0 lambda B = 7*pi; % assign parameters

Bo = 0.3*B;

for m = 0:1:200 for n = 0:1:200

if(((m-100)*(m-100)+(n-100)*(n-100))<1600)

Q1(m,n)=exp(j*2*pi*(Bo*((((m-100)/40)^2+((n-100)/40)^2)^3)+B*((((m-100)/40)^2…

+((n-100)/40)^2)^3)*log((((m-100)/40)^2+((n-100)/40)^2)^0.5 +1e-8)));

Q2(m,n)=exp(j*2*pi*(Bo*((((m-100)/40)^2+((n-100)/40)^2)^3)+B*((((m-100)/40)^2…

+((n-100)/40)^2)^3)*log((((m-100)/40)^2+((n-100)/40)^2)^0.5 +1e-8)));

Q3(m,n)=exp(j*2*pi*(Bo*((((m-100)/40)^2+((n-100)/40)^2)^3)+B*((((m-100)/40)^2…

+((n-100)/40)^2)^3)*log((((m-100)/40)^2+((n-100)/40)^2)^0.5 +1e-8)));

Q4(m,n)=exp(j*2*pi*(Bo*((((m-100)/40)^2+((n-100)/40)^2)^3)+B*((((m-100)/40)^2…

+((n-100)/40)^2)^3)*log((((m-100)/40)^2+((n-100)/40)^2)^0.5 +1e-8)));

end end

end

for m = 0:1:200

for n = 0:1:200

if(((m-100)*(m-100)+(n-100)*(n-100))<1600)

Q1(m,n)=Q1(m,n)*exp(j*2*pi*((((m-100)/40)^2+((n-100)/40)^2)^3)*0);

Q2(m,n)=Q2(m,n)*exp(j*2*pi*((((m-100)/40)^2+((n-100)/40)^2)^3)*0.5);

Q3(m,n)=Q3(m,n)*exp(j*2*pi*((((m-100)/40)^2+((n-100)/40)^2)^3)*1);

Q4(m,n)=Q4(m,n)*exp(j*2*pi*((((m-100)/40)^2+((n-100)/40)^2)^3)*2);

end end

end

otf1=conv2(Q1,conj(Q1));

otf2=conv2(Q2,conj(Q2));

otf3=conv2(Q3,conj(Q3));

otf4=conv2(Q4,conj(Q4));

for n = 1:80

mtf1(n)=abs(otf1((n+198),199));

mtf2(n)=abs(otf2((n+198),199));

mtf3(n)=abs(otf3((n+198),199));

mtf4(n)=abs(otf4((n+198),199));

end

% Normalization mtf1= mtf1/max(mtf1);

mtf2= mtf2/max(mtf2);

mtf3= mtf3/max(mtf3);

mtf4= mtf4/max(mtf4);

x1=0:(2/80):(2-2/80);

plot(x1,mtf1,x1,mtf2,x1,mtf3,x1,mtf4)

Appendix C: MATLAB Source Codes (CONTINUED) C.4: Maximal on-axis intensity versus BB46

.

% --- TOPIC: Maximal on-axis intensity versus BB46 --- %

% --- Drafted by Chih-Yun Chan, 1st version, test ok --- %

% --- Used in Section 2 --- %

u1 = zeros(161,1); % Field I1 = zeros(161,1); % Intensity

maxi = zeros(81,1); % Max. intensity v.s. B46

for B46 = 0: 1 : 80 for B26 = 0 : 1 : 160 for x = 0 : 0.001 : 1

y1 = x*exp(j*2*pi*1*((B26./20 -4)*(x^2)+(B46./20-3.5)*(x^4)+(x^6)));

u1(B26+1)= u1(B26+1)+ 0.001*y1;

end end

I1 = u1.*conj(u1)*4*pi*pi/10;

maxi(B46+1)=max(I1);

u1 = zeros(161,1);

I1 = zeros(161,1);

end

k= -3.5:0.05:0.5;

plot(k,maxi)

Appendix C: MATLAB Source Codes (CONTINUED)

C.5: Plot of the Strehl ratio versus B26 for different aberration scaling factors f, with B46 = -1.5

% --- TOPIC: Strehl ratio versus B26 for different f, with BB46 = -1.5 --- %

% --- Drafted by Chih-Yun Chan, 1st version, test ok --- %

% --- Used in Section 2 --- %

u0 = zeros(101,1); u1 = zeros(101,1); u2 = zeros(101,1); u3 = zeros(101,1);

u4 = zeros(101,1); u5 = zeros(101,1); u6 = zeros(101,1); % field with ideal lens

u02 = zeros(101,1); u12 = zeros(101,1); u22 = zeros(101,1); u32 = zeros(101,1);

u42 = zeros(101,1); u52 = zeros(101,1); u62 = zeros(101,1); % field with the filter

I0 = zeros(101,1); I1 = zeros(101,1); I2 = zeros(101,1); I3 = zeros(101,1);

I4 = zeros(101,1); I5 = zeros(101,1); I6 = zeros(101,1); % intensity with ideal lens

I02 = zeros(101,1); I12 = zeros(101,1); I22 = zeros(101,1); I32 = zeros(101,1);

I42 = zeros(101,1); I52 = zeros(101,1); I62 = zeros(101,1); % intensity with the filter

B46=-1.5; % assign parameters B = 3.9*pi;

Bo = 0.55*B;

for B26 = 0 : 1 : 100 % doing the integration for x = 0 : 0.001 : 1

y0 = x*exp(j*2*pi*0.5*((B26./100)*(x^2)+B46*(x^4)+(x^6)));

y02 = x*exp(j*2*pi*0.5*((B26./100)*(x^2)+B46*(x^4)...

+(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u0(B26+1)= u0(B26+1)+ 0.001*y0;

u02(B26+1)= u02(B26+1)+ 0.001*y02;

y1 = x*exp(j*2*pi*1*((B26./100)*(x^2)+B46*(x^4)+(x^6)));

y12 = x*exp(j*2*pi*1*((B26./100)*(x^2)+B46*(x^4)...

+(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u1(B26+1)= u1(B26+1)+ 0.001*y1;

u12(B26+1)= u12(B26+1)+ 0.001*y12;

y2 = x*exp(j*2*pi*2*((B26./100)*(x^2)+B46*(x^4)+(x^6)));

y22 = x*exp(j*2*pi*2*((B26./100)*(x^2)+B46*(x^4)...

+(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u2(B26+1)= u2(B26+1)+ 0.001*y2;

u22(B26+1)= u22(B26+1)+ 0.001*y22;

y3 = x*exp(j*2*pi*3*((B26./100)*(x^2)+B46*(x^4)+(x^6)));

y32 = x*exp(j*2*pi*3*((B26./100)*(x^2)+B46*(x^4)...

+(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u3(B26+1)= u3(B26+1)+ 0.001*y3;

u32(B26+1)= u32(B26+1)+ 0.001*y32;

y4 = x*exp(j*2*pi*4*((B26./100)*(x^2)+B46*(x^4)+(x^6)));

y42 = x*exp(j*2*pi*4*((B26./100)*(x^2)+B46*(x^4)...

+(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u4(B26+1)= u4(B26+1)+ 0.001*y4;

u42(B26+1)= u42(B26+1)+ 0.001*y42;

y5 = x*exp(j*2*pi*5*((B26./100)*(x^2)+B46*(x^4)+(x^6)));

y52 = x*exp(j*2*pi*5*((B26./100)*(x^2)+B46*(x^4)...

+(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u5(B26+1)= u5(B26+1)+ 0.001*y5;

u52(B26+1)= u52(B26+1)+ 0.001*y52;

y6 = x*exp(j*2*pi*6*((B26./100)*(x^2)+B46*(x^4)+(x^6)));

y62 = x*exp(j*2*pi*6*((B26./100)*(x^2)+B46*(x^4)...

+(x^6)))*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

u6(B26+1)= u6(B26+1)+ 0.001*y6;

u62(B26+1)= u62(B26+1)+ 0.001*y62;

end end

I0 = u0.*conj(u0)*4*pi*pi/10;

I1 = u1.*conj(u1)*4*pi*pi/10;

I2 = u2.*conj(u2)*4*pi*pi/10;

I3 = u3.*conj(u3)*4*pi*pi/10;

I4 = u4.*conj(u4)*4*pi*pi/10;

I5 = u5.*conj(u5)*4*pi*pi/10;

I6 = u6.*conj(u6)*4*pi*pi/10;

I02 = u02.*conj(u02)*4*pi*pi/10;

I12 = u12.*conj(u12)*4*pi*pi/10;

I22 = u22.*conj(u22)*4*pi*pi/10;

I32 = u32.*conj(u32)*4*pi*pi/10;

I42 = u42.*conj(u42)*4*pi*pi/10;

I52 = u52.*conj(u52)*4*pi*pi/10;

I62 = u62.*conj(u62)*4*pi*pi/10;

k = 0 : 0.01 : 1.0;

plot(k,I12,k,I22,k,I32,k,I42,k,I52,k,I62)

Appendix C: MATLAB Source Codes (CONTINUED) C.6: The Strehl ratio as a function of aberration scaling factor f, with

BB26

=0.6 and B

46

= -1.5.

% --- TOPIC:Strehl ratio as a function of f, with B26 =0.6 and B46 = -1.5. --- %

% --- Drafted by Chih-Yun Chan, 2nd version, test ok --- %

% --- Used in Section 2 --- %

u1 = zeros(101,1); % field with ideal lens u2 = zeros(101,1); % field with the phase filter I1 = zeros(101,1); % intensity with ideal lens I2 = zeros(101,1); % intensity with the phase filter

B46=-1.5; % assign parameters B = 5*pi;

Bo = 0.3*B;

for f = 0 : 1 :100 % doing the integration for x = 0 : 0.001 : 1

y1 = x*exp(j*2*pi*f/10*(0.6*(x^2)+B46*(x^4)+(x^6)))…

*[exp(j*2*pi*(Bo*(x^6)+B*(x^6)*log(x+1e-6)))];

y2 = x*exp(j*2*pi*f/10*(0.6*(x^2)+B46*(x^4)+(x^6)));

u1(f+1)= u1(f+1)+ 0.001*y1;

u2(f+1)= u2(f+1)+ 0.001*y2;

end end

I1 = u1.*conj(u1)*4*pi*pi/10;

I1= I1/max(I1); % Normalization I2 = u2.*conj(u2)*4*pi*pi/10;

I2= I2/max(I2); % Normalization

k = 0 : 0.1 : 10;

plot(k,I1,k,I2)

Appendix C: MATLAB Source Codes (CONTINUED)

C.7: 1) Relation among the transverse gain and the radius of the first zone.

2) Relation among the Strehl ratio and the radius of the first zone

-- For CD

% --- TOPIC: SR v.s. radius a; GT v.s. radius a (for CD) --- %

% --- Drafted by Chih-Yun Chan, 2nd version, test ok --- %

% --- Used in Section 4 --- %

double NUM; % number of points in the interval a: 0.4~0.6 double a(NUM+1); % radius of the 1st zone

double b(NUM+1); % radius of the 2nd zone double t1; % transmission of zone 1 double t2; % transmission of zone 3 double phi; % phase difference of the zones double I0(NUM+1); double a0(NUM+1); double b0(NUM+1);

double I1(NUM+1); double a1(NUM+1); double b1(NUM+1);

double I2(NUM+1); double a2(NUM+1); double b2(NUM+1);

double uF(NUM+1);

double SR(NUM+1);

double GT(NUM+1);

double GA(NUM+1);

t1= 1; % assign the parameters t2= 1;

phi=0.4 NUM = 200;

for i = 0 : NUM

a(i+1) = 0.4 + (0.2/NUM).*i;

b(i+1) = sqrt(1-a(i+1).*a(i+1));

I0(i+1) = exp(j*0*pi)*(t1.*(a(i+1).^2)) + t2.*exp(j*phi*pi)*(1-(b(i+1).^2));

a0(i+1) = real(I0(i+1));

b0(i+1) = imag(I0(i+1));

I1(i+1) = 0.5*exp(j*0*pi)*(t1.*(a(i+1).^4)) + t2.*0.5*exp(j*phi*pi)*(- (b(i+1).^4) +1);

a1(i+1) = real(I1(i+1));

b1(i+1) = imag(I1(i+1));

I2(i+1) = (1/3)*exp(j*0*pi)*(t1.*(a(i+1).^6))+ t2.*(1/3)*exp(j*phi*pi)*(- (b(i+1).^6) +1);

a2(i+1) = real(I2(i+1));

b2(i+1) = imag(I2(i+1));

% the displacement of the focus

uF(i+1) = 2*(a1(i+1).*b0(i+1)-a0(i+1).*b1(i+1))./( (a2(i+1).*a0(i+1)+b2(i+1).*b0(i+1))…

- (a1(i+1).*a1(i+1)+b1(i+1).*b1(i+1)) );

% the Strehl ratio

SR(i+1) = a0(i+1).*a0(i+1) + b0(i+1).*b0(i+1)…

- uF(i+1).*(a0(i+1).*b1(i+1)-a1(i+1).*b0(i+1));

% the transverse gain

GT(i+1) = 2*( (a1(i+1).*a0(i+1)+b0(i+1).*b1(i+1))…

- uF(i+1).*(-a2(i+1).*b0(i+1)+a0(i+1).*b2(i+1)) )./SR(i+1);

% the axial gain

GA(i+1) = 12*( (a2(i+1).*a0(i+1)+b0(i+1).*b2(i+1))…

- (a1(i+1).*a1(i+1)+b1(i+1).*b1(i+1)) )./SR(i+1);

end

m = 0.4: (0.2/NUM): 0.6;

subplot(2,1,1); plot(m,uF);

subplot(2,1,2); plot(m,SR);

Appendix C: MATLAB Source Codes (CONTINUED)

C.8: 1) Relation among the transverse gain and the radius of the first zone.

2) Relation among the Strehl ratio and the radius of the first zone

-- For DVF

% --- TOPIC: SR v.s. radius a; GT v.s. radius a (for DVD) --- %

% --- Drafted by Chih-Yun Chan, 2nd version, test ok --- %

% --- Used in Section 4 --- %

double NUM; % number of points in the interval a: 0.4~0.6 double a(NUM+1); % radius of the 1st zone

double b(NUM+1); % radius of the 2nd zone double t1; % transmission of zone 2 double t2; % transmission of zone 4 double phi; % phase difference of the zones double I0(NUM+1); double a0(NUM+1); double b0(NUM+1);

double I1(NUM+1); double a1(NUM+1); double b1(NUM+1);

double I2(NUM+1); double a2(NUM+1); double b2(NUM+1);

double uF(NUM+1);

double SR(NUM+1);

double GT(NUM+1);

double GA(NUM+1);

c= 0.7; % assign parameters t1=1;

t2=1;

phi=0.35;

NUM = 200;

for i = 0 : NUM

a(i+1) = 0.4 + (0.2/NUM).*i;

b(i+1) = sqrt(1-a(i+1).*a(i+1));

I0(i+1) = t1.*exp(j*0*pi)*((c.*b(i+1)).^2-(c.*a(i+1)).^2) + t2.*exp(j*phi*pi)*(1-(c).^2);

a0(i+1) = real(I0(i+1));

b0(i+1) = imag(I0(i+1));

I1(i+1) = (1/2)*(t1.*exp(j*0*pi)*((c.*b(i+1)).^4-(c.*a(i+1)).^4) + t2.*exp(j*phi*pi)*(1-(c).^4));

a1(i+1) = real(I1(i+1));

b1(i+1) = imag(I1(i+1));

I2(i+1) = (1/3)*(t1.*exp(j*0*pi)*((c.*b(i+1)).^6-(c.*a(i+1)).^6) + t2.*exp(j*phi*pi)*(1-(c).^6));

a2(i+1) = real(I2(i+1));

b2(i+1) = imag(I2(i+1));

% the displacement of the focus

uF(i+1) = 2*(a1(i+1).*b0(i+1)-a0(i+1).*b1(i+1))…

./( (a2(i+1).*a0(i+1)+b2(i+1).*b0(i+1)) - (a1(i+1).*a1(i+1)+b1(i+1).*b1(i+1)) );

% the Strehl ratio

SR(i+1) = a0(i+1).*a0(i+1) + b0(i+1).*b0(i+1) - uF(i+1).*(a0(i+1).*b1(i+1)-a1(i+1).*b0(i+1));

% the transverse gain

GT(i+1) = 2*( (a1(i+1).*a0(i+1)+b0(i+1).*b1(i+1))…

- uF(i+1).*(-a2(i+1).*b0(i+1)+a0(i+1).*b2(i+1)) )./SR(i+1);

% the axial gain

GA(i+1) = 12*( (a2(i+1).*a0(i+1)+b0(i+1).*b2(i+1)) … -(a1(i+1).*a1(i+1)+b1(i+1).*b1(i+1)) )./SR(i+1);

end

m = 0.4: (0.2/NUM): 0.6;

%subplot(2,1,1);

plot(m,uF);

%subplot(2,1,2); plot(m,SR);

Appendix C: MATLAB Source Codes (CONTINUED) C.9: The transverse intensity distributions -- For CD

% --- TOPIC: the transverse intensity distributions (for CD) --- %

% --- Drafted by Chih-Yun Chan, 2nd version, test ok --- %

% --- Used in Section 4 --- %

v = (0:0.01:10)'; % radius distance double a(3); % radius of the 1st zone double b(3); % radius of the 2nd zone double U(1001,3); % field with the clear aperture double Ut(1001,3); % field with the designed filter double I(1001,3); % intensity with the clear aperture double It(1001,3); % intensity with the designed filter

a(1)=0.5; a(2)=0.57; a(3)=0.6;

b(1)=sqrt(1-a(1).*a(1)); b(2)=sqrt(1-a(2).*a(2)); b(3)=sqrt(1-a(3).*a(3));

warning off;

for k = 1:3 for i = 0 : 1000

U(i+1,k) = (2./v(i+1))*besselj(1,v(i+1));

Ut(i+1,k) = (2./v(i+1))*( a(k).*besselj(1,(a(k).*v(i+1)))…

+ exp(j*0.1*pi).*( besselj(1,v(i+1)) - b(k).*besselj(1,(b(k).*v(i+1))) ) );

end end

for k= 1:3

I(1:1001,k) = U(1:1001,k).*conj(U(1:1001,k));

It(1:1001,k) = Ut(1:1001,k).*conj(Ut(1:1001,k));

end

subplot(2,1,1); plot(v,I(:,1),v,I(:,2),v,I(:,3));

subplot(2,1,2); plot(v,It(:,1),v,It(:,2),v,It(:,3));

warning on;

Appendix C: MATLAB Source Codes (CONTINUED) C.10: The transverse intensity distributions -- For DVD

% --- TOPIC: the transverse intensity distributions (for DVD) --- %

% --- Drafted by Chih-Yun Chan, 2nd version, test ok --- %

% --- Used in Section 4 --- %

v = (0:0.01:10)'; % radius distance double a(3); % radius of the 1st zone double b(3); % radius of the 2nd zone double U(1001,3); % field with the clear aperture double Ut(1001,3); % field with the designed filter double I(1001,3); % intensity with the clear aperture double It(1001,3); % intensity with the designed filter

a(1)=0.5; a(2)=0.57; a(3)=0.6;

b(1)=sqrt(1-a(1).*a(1)); b(2)=sqrt(1-a(2).*a(2)); b(3)=sqrt(1-a(3).*a(3));

c=0.7;

warning off;

warning off;