More Efficient use of Recursion Formulas to create arrays and to flatten arrays into data vectors

[CS_10222121]

Hi all:

I am looking for help/tips regarding more efficient/practical ways in MathCAD 15 to perform recursive calculations as it relates to creation of arrays, vectorizing the arrays (row by row), and then augmenting these data vectors into a final array.  I have attached MathCAD 15 file to illustrate my questions.

The attached file is intended to perform a least-squares fitting of optical wavefront data (optical path difference data) to a set of Zernike polynomials.  I have tried to describe each portion of the script so that it is easier to understand what I am trying to accomplish.

The first portion of the script creates two arrays (rho and theta) that contain the respective values over the unit circle.

As shown, the next section explicitly lists the mathematical form for of the first 36 Zernike polynomials and evaluates them at each coordinate position.  Once evaluated, individual arrays are created for each, evaluated Zernike polynomial so that in the end we have 36 individual arrays of the proper dimensions.  As many of you are probably aware, the Zernike polynomials are an orthogonal set of functions that can be generated up to any order using a recursion formula.  You will note that I have typed out the radial portion of the Zernike polynomial and the conditional rules for angular portion.  For this portion of the routine what I am hoping is if someone can help me create a solve block or iterative loop that allows me to only write out the recursion formulas for the Zernike polynomials if I provide a range for the n and integer values.

The portion of the routine takes the 36 individual Zernike arrays, converts them into data vectors (row by row), and then augments them into a single array.  The script in this section borrows from a method that I saw somebody show on the PTC community board a couple of years ago.  This approach does work, but again was wondering if there was a way to to perform these operations more efficiently.

The next section simply reads in the data file and then flattens it into a data vector.  The last section performs a series of matrix operations to make the arrays the correct dimensions and to calculate the coefficients for each individual Zernike term.

In brief, I am wondering if someone can help me more efficiently write out the second and third sections of this script by using the recursion formulates for Zernike polynomials and a more efficient way to flatten the arrays into data vectors.  Thanks in advance.

[terryhendicott]

Hi,

Taken a while to answer but the enclosed worksheet will help.

Cheers

Terry

[CS_10222121]

Thanks Terry.  I appreciate the response and work.  I will review and let you know if I have any questions.

[terryhendicott]

Hi,

Have made a correction and some additional information.  See items in yellow.

If you set the 3D plot to the following setting you get a plan view of disk in exact orientation as the Wikipedia disks.

Cheers

Terry

[terryhendicott]

Hi,

I have a question regarding the 517 x 517 pixel square and Zernike polynomials defined on the unit circle.  Before I blast away and show you how to file the 517x517 pixel grid.

Zernike Polynomial radial function is only applicable to radius "rho" from 0.0 to 1.0. The unit circle. In the case of a square 517 x 517 grid your calculations have a "rho" in the corners and adjacent areas of up to 1.414 being the diagonal (radius) of the unit square. Do you want to do this or do you want a conformal mapping square to unit circle with "rho" <= 1.0? Below describes one such mapping.   Doing some research to get a reversable method of mapping square to circle.

Have plotted them so you see what I mean.

Cheers

Terry

[CS_10222121]

Hi Terry

thanks again for the reply and interest in this problem.  The reason for the 517x517 pixel square in my example script is because the wavefront data I am trying to fit had these dimensions.  The data itself is over a circular aperture that can conform to a unit circle.  You are correct in that the Zernike polynomial radial function as defined is only applicable to a radius from 0 to 1..  As you will note from my example script, when I wrote out explicitly the first 36 Zernike polynomials and evaluated them at each value of rho and theta, I did define a rmax value that was used to normalize and evaluate each Zernike polynomial at the maximum r value.  I effectively defined a max radius and then only evaluated each zernike polynomial at each rho and theta value if it was less than or equal to the max radius.  While I haven't delved into the details for the script your sent, the resultant images look like you "extracted" a unit circle from the initial square array corresponding to the individual rho and theta values over the unit circle.  Is this correct?  If so, this seems like a good idea.  let me know if I am misunderstanding your question.  Thanks. 

[terryhendicott]

Hi,

My bad.  Did not look at your file closely enough.  Understand now you "extract" a circle in the middle of the 517x517 square and set 0.0 as value outside this middle circle.

The mapping I asked about "squashes and deforms" the square into a unit circle.  This is clearly not what you want done.  I have programmed the reversal method for "unsquashing and undeforming" from circle to square but that is mute now.

Will think about the problem again.

Cheers

Terry

[terryhendicott]

Hi,

Done!

Now coordinates are by pixel and values outside the focus are zero.

Warning the resulting "file" matrix is large.  I have 128 GB of ram so can have really large matrices.

If you have problems reduce the maximum order of polynomials

Cheers

Terry

[ttokoro]

Prime only plot 40000 points for each color of plot points.

[Werner_E]

---

Warning the resulting "file" matrix is large.  I have 128 GB of ram so can have really large matrices.

We are using a 32-bit application here, which by design can theoretically only use up to 4GB of RAM. 😉

But a fast processor sure helps cutting the calculation time.

[CS_10222121]

Hi Terry

Thanks for the new script!  I ran the script and, from what I can determine, it does evaluate the user selected # of Zernike polynomials and then "flattens" each evaluated Zernike polynomial into a data vector and joins each of these data vectors to create a single array.  It is then straight forward to perform some matrix algebra with the data to perform the fitting.  One additional thing I forgot to mention and appears in the original file I sent for each individual Zernike polynomial is a normalization constant.  The definition of the Normalization constant for Zernike polynomials is

Nn,m = sqrt((2n+2)/(1+delta)), where delta = 1 for m = 0 and delta = 0 for m not equal to zero.

In looking at your script it would appear that I could define this normalization constant after the definition of the radial polynomial and then combine this with your definition of the full Zernike polynomial. Not sure how to enter the logic You might have a better suggestion.  thanks.

[terryhendicott · Accepted solution]

Hi,

Have included a worksheet that shows the difference between:-

a) calculating by formulas Zj(rho,theta) that needs a Normalization constant, and

b) calculating using the direct polynomial definition Znm(rho,theta) that does not need the Normalization constant. 

The Normalization worksheet shows the difference between a) and b)

I have made some alterations to the original file that stops recreation of the "file" matrix when you change the plot above.

Cheers

Terry

[CS_10222121]

Thanks Terry.  Script runs well.  Appreciate all of your assistance with this matter.

[AlanStevens]

Not that I know what I'm doing here (unlike Terry), but this looked like it could be fun! So, perhaps the attached might be of some use.

I get the following images, which are similar to those on Wikipedia, though rotated:

[terryhendicott]

Hi Alan,

Do not profess to know what I am doing.  Like you do it for fun.

Cheers

Terry