Archive for March, 2007

Papers

Tuesday, March 27th, 2007

I've started using Papers to complement EndNote for keeping track of research papers that I've read. Think iPhoto for papers. In the past I've manually named papers that I've downloaded and tried to file them away. Except folders aren't very good metadata containers, because a paper will often need multiple tags. Papers handles all of the metadata for you and stores the files away on your disk. It has a built-in fulls screen view and uses PubMed to lookup random PDF files that you give it.

Why not just use EndNote? Its interface, file management and metadata capabilities are like a blast from the past circa MacOS 9. The programmers have probably been too busy rolling around in their money to actually make it decent.

Papers still needs some work though. It can't export to Word, which is why I still use both Papers and EndNote. It would be nice if it automatically pulled the paper title and looked it up in PubMed for you for identification. Speaking of which, they should hook it up with something like Google Scholar that is used by everyone. PubMed is pretty broad but definitely isn't used by the math and physics people. So Papers shows some promise and I'm using it for about 75% of my needs right now so stay tuned.

Beer Brewing

Sunday, March 18th, 2007

I've been brewing beer on a very irregular basis since the summer of 2002 and have made 7 batches thus far. Everything that I've learned has been based upon a few used books, friends and family and keeping meticulous records.

Why Brew Beer?

Jordan and I started brewing beer for two reasons. First, we had just graduated from high school and were not of legal drinking age yet. Second, good beer is expensive. A 12-pack of Sierra Nevada will set you back about the same amount in dollars and brewing seemed like a reasonable way to make any flavor of beer for the price of Pabst Blue Ribbon.

A more recent motivation to learn more about the basic processes involved is that (Saccharomyces cerevisiae) is the probably most intensely studied eucaryotic cell in history. It's genome composed of about 6000 genes was completely mapped in 1997 and there is plenty of information out there. There's also plenty of information about the yeast life cycle and ideal growth conditions. Brewing beer is just a practical application of insanely complicated biochemistry to create a tasty beverage and reduce your chance of scurvy.

Fermentation Temperature Control

Bottling

Bottling has a much lower initial cost and time investment than kegging, however it will quickly become apparent that the investment in a good kegging system will pay for itself very quickly. Between the 2-4 hours required to setup and bottle 5 gallons of beer and the ongoing cost of replacement bottles for the ones that grew all kinds of interesting things, kegging is cheaper in the long run.

With that said, you will need the following to bottle one 5-gallon batch of homebrew:

  • (60) Relatively clean 12 oz. glass bottles (no screw-on caps)
  • (60) Bottle caps
  • (1) Hand-held bottle capper
  • (1) Beer filled glass carboy
  • (8 feet) 1/4" ID clear hose
  • (1) Bottling Tube
  • (1) Bottle scrubbing brush
  • (1 cup) Corn sugar
  1. Make the priming solution
    • Bring the corn sugar and 16 ounces of water to a boil. Stir until the sugar is dissolved. Cover, remove from heat and let cool until needed again.
  2. Sanitize the bottles, caps, bottling tube and clear hose
    • Fill a sink with cold water and iodine cleaning solution per the instructions on the bottle
    • Quickly scrub the insides of the bottles with your scrubbing brush to remove any grime
    • Place them upside down on paper towels and allow to dry as long as possible.
  3. Prepare for bottling
    • Place the carboy on a very high surface to generate enough surface for racking to the bottles.
    • Add the priming solution to the carboy. It will diffuse fairly uniformly by the time that you fill the bottles.
    • Place the clear tube in the carboy and lower to within 1-2" of the yeast sediment line
    • Attach the bottling tube to the clear hose, hold it above the carboy. Insert the valve side into your mouth and press with your tongue. Suck until enough beer is in the tube to start a siphon
    • Gather as many bottles as is convenient on the ground for easy access. Sit with them.
  4. Bottle!
    • Press the bottling tube to the bottom of a bottle and fill until the beer is at the very top of the bottle. The displacement of the bottling tube will correct the apparent fill level (Be sure to fill the bottles with enough beer, or else there will be excess oxygen in the bottles that will allow hearty aerobic bacteria to grow quickly and take over your beer more quickly than normal)
    • Bottle until all of the bottles are full or you run out of beer. Don't be afraid to leave a small amount of good beer in the carboy to minimize the amount of yeast that ends up in each bottle.
  5. Cap!
    • Carefully place the bottle caps on all of your bottles and then cap them. Once you get the technique down it should go fairly quickly.

Kegs

Keeping it cold: Building a kegerator on the cheap.

to be continued...

Write your own Wedding Script

Sunday, March 18th, 2007

Our friend Andrew performed the ceremony at our wedding. We wrote up the script together beforehand using Writely. We wanted something to the point and the ceremony lasted about 7 minutes I think. Without further ado, here's the script:

Introduction

Minister: Dearly beloved, we are gathered here together today to join together this man and this woman in holy matrimony. We are here to celebrate their union and to honor their commitment to one another. Today, J and C proclaim their love to the world and we rejoice with them.

In marriage, we give ourselves freely and generously into the hands of the one we love, and in doing so, each of us receives the love and trust of the other as our most precious gift. But even as that gift is shared by two people who are in love, it also touches the friends and family members who in various ways support and contribute to the relationship. All of you are J's and C’s community, and each of you has played some part in bringing them to this moment. This is why gathering as a community is such an important part of a wedding ceremony. J and C are now taking a new form as a married couple, and in this form, they become part of our community in a new way.

I'm proud to say that i am part of J and C's community, and have had the privilege of witnessing these two interact over the last three years. From the beginning it was evident that their relationship was somehow different than those of their peers; something about how they looked at each other, or the care and compassion that they showed for each other from the very beginning.

I recall one particular incident when the three of us were driving in a car together. It was after their second year in college and summer was just beginning. C was moving to her new apartment, and we were helping her move. I asked, why aren't you living together? It seemed like an obvious question. They both looked at each other and didn't have an answer for me. It seemed like they just hadn't thought about it. At that point, it was already obvious to me and probably many of you that they were meant for each other. And sure enough, they moved in together one year later. the following year they went to Yosemite for C's birthday and returned newly engaged. Hearing of their engagment was much the same as hearing that they had moved in together, it just seemed like it was inevitable

A vast, unknown future stretches out before you. The future, with its hopes and disappointments, its joys and its sorrows, is hidden from your eyes. But it is a great tribute to your belief in each other that you are willing to face those uncertainties together. May the pure, simple love with which you join hearts and hands today never fail, but grow deeper and surer with every year you spend together.

The Marriage Vows

J and C, we are here to remember and rejoice with you and to recount with one another that it is love that guides us on our path, and to celebrate as you begin this journey together. It is in this spirit that you have come here to today to exchange these vows.

Minister: (To the groom) J, do you take this woman to be thy wedded wife? To have and to hold from this day forward, for better or for worse, for richer, for poorer, in sickness and in health, to love and to cherish; from this day on so long as you both shall live?

Groom: I do.

Minister: (To the bride) C, do you take this man to be thy wedded husband? To have and to hold from this day forward, for better or for worse, for richer, for poorer, in sickness and in health, to love and to cherish; from this day on so long as you both shall live?

Bride: I do.

The Rings

Minister: The wedding ring is an unbroken circle symbolizing unending and everlasting love.

J, have you a token of your love for C?

(A hands Joey the ring)

Minister: C, have you a token of your love for J?

(M hands C the ring)

Minister: Traditionally, the passage to the status of husband and wife is marked by the exchange of rings. These rings are a symbol of the unbroken circle of love. Love freely given has no beginning and no end. Love freely given has no giver and no receiver - for each is the giver and each is the receiver. May these rings remind you always of the vows you have taken here today.

Minister: J place this ring on C's finger and repeat after me:

This ring, a gift for you, symbolizes my desire that you be my wife from this day forward.

Groom: (repeats)

Minister: C place this ring on J's finger and repeat after me:
This ring, a gift for you, symbolizes my wish that you be my husband from this day forward.

Bride: (repeats)

Ending

Minister: J and C, you have given and pledged your promises to each other, and have declared your everlasting love by giving and receiving rings. By the authority vested in me by the State of California, I now pronounce you husband and wife.

(To the groom) You may kiss the bride.

(cheering)

Minister: Ladies and Gentlemen, I present to you .

Cats

Sunday, March 18th, 2007

Flea Medicine

The flea medicine business is kind've funny. Medicine for your 10 pound cat costs as much as for your 150 pound wildabeast even though the only difference is the volume of the dose. Here's the info on Frontline plus for the dog medicine.

Pet Weight (lbs) Recommended Dose (ml) Price for 6 doses
0-22 0.67 62.99
23-44 1.34 67.95
45-88 2.68 63.95
89-132 4.52 73.95

So what we do is buy the big dog dose, transfer it to a small glass bottle and store it in the fridge. When it's time to give the cats their medicine we use a syringe (no needle!) to dispense the medicine. For two cats at 0.5 ml per dose, 4.5 ml will last 5 months, so you're talking about $1 per dose rather than $10 per dose.

Costs/Reviews

Our cat sherlock had an endoscopy to remove a needle that he swallowed. It was performed in Capitola, CA at Pacific Veterinary Specialists Emergency Services on a Saturday night and cost about $1000.

In July 2007 Sherlock started having trouble urinating and we saw that he had blood in his urine. We took him to Adobe Pet Hospital because they were nearby and had drop-in appointments. They had to insert a catheter and kept him for two nights. The overall cost was about $1300. Overall it was pretty poor service and it didn't seem like the veterinarians or staff cared about how our cat did. It was in stark contrast with our previous encounter at PVES (see above). I wouldn't go there again.

Later in July, Sherlock had serious trouble going to the bathroom again. They weren't able to get the catheter in because it was so clogged, so he ended up going into surgery (perineal urethrostomy) to open up his urethra. The surgery went great and he has been doing really well since. It was expensive, but it worked out to be close to the same cost of the prior hospital stays when he had to have a catheter (~$3000) and the fact that he can go to the bathroom is infinitely worth it.

Doggy Daycare

We take Lola (our dog) to Klub K9 Pet Center in Sunnyvale and highly recommend it. It's a little over $20/day if you buy a 20-pack, otherwise it's about $40 per day. Their staff is great, they have four huge play areas for the dogs, and they do cage-free boarding.

Backpacking

Sunday, March 18th, 2007

Here's a backpacking list that I originally got from my dad.

Clothing
 Kitchen Personal Gear
2 pairs of socks stove + fuel sunscreen
1 pair underwear pots + utensils deet
1long sleeve shirt matches + lighter camera + batteries
1 pair of long pants + zipoffs swiss army knife pencil + paper
warm outerwear lexan cup books
rain gear folding water bucket precscriptions
warm hat water filter + tablets survival kit*
water shoes scouring pad first aid kit **
  bowl + spoon head light + batteries
wine repair kit ***
bleach small daypack / fannypack
hygiene sleeping trailhead hike in
toothbrush +paste + floss sleeping bag hiking poles
toilet paper + shovel sleeping pad sunglasses
wash cloth p bottle water bottle + drink mix
camp suds tent/ground clot sun hat
lip balm
sunscreen
bandanna
map
id + credit card + cash
 survival kit
medical kit
repair kit
cup steristrips bobby pins
waterproof matches ace bandage super glue
firestarter moleskin spare bulbs
whistle a&d ointment duct tape
space blanket ibuprofen eyeglass repair kit
2 yards toilet paper benadryl sewing kit
50 feet nylon cord immodium
note paper + pencil alcohol wipes
tiny knife tweezer
tiny light 4x4 gauze pads
3 bouillon cubs bandaids
3 teabags bacitracin
3 hard candies second skin

Travel

Sunday, March 18th, 2007

See Camping and Backpacking for packing checklists and other information.

Conquered

  • San Diego and Healdsburg
  • Costa Rica
  • Yosemite National Park, CA
  • Switzerland

Unexplored

Near Home

*Hawaii
*Yellowstone National Park, Wyoming
*Mount Rushmore, South Dakota
*Niagara Falls, ON, Canada
*New York, New York
*Boston, Massachusetts
*The Arch, St Louis, Missouri
*New Orleans, Louisiana
*Mammoth Cave, Kentucky

Proposed Trips

*Southwestern Extravaganza
**Yosemite, CA
**King's Canyon
**Bryce
**Zion
**Las Vegas, Nevada
**Hoover Dam, Boulder City, Nevada
**Grand Canyon, AZ
**Meteor Crater, AZ
**Petrified Forest and Painted Desert, AZ
*Honneymoon in the Bahamas

An RSVP Application in Ruby on Rails

Friday, March 9th, 2007

I wrote an RSVP program for our wedding in Ruby on Rails. It includes live search and editing using AJAX which was ridiculously easy with Rails. It also has a shadow table to which all changes are written in case one of your guests gets all drunk and uppity while signing up. Here is the source code, including the pretty graphics. Feel free to use it for anything and let me know if you find it useful.

Schema

The database schema is included. You can import it with:

mysql -u username dbname < schema.mysql

RSVP rails code and schema

Create a Nice Looking Skin for Mediawiki

Saturday, March 3rd, 2007

The default installation of MediaWiki is extremely functional but not too easy on the eyes. There are some nice examples of skins online. Particularly good looking examples are the Hula and Beagle projects. The Beagle guys are really nice and released their skin, which I used as an example for this website's previous incarnation. There were a couple of bugs on the Beagle skin that I fixed and simplified for this site.

In order to customize the skin you will need to walk through the code. Here's a brief description of the files.

  • Beagle.php: Specifies where everything fits together. If you want to move the navigation bar, eliminate the searchbar or enable Google Analytics this is the place to do it.
  • main.css: Specifies everything from text style, size and color to the look of the navigation bar.
  • logo.png: The Dogully tree logo up at the top of your screen
  • title_bg.png: The orange gradient that sits behind the logo. Specified in main.css

Calculate Damping Coefficients in Matlab

Friday, March 2nd, 2007

When working with experimental data you'll often run into damped oscillations. Most of the time you can eyeball the data and see that the disturbance is gone before it affects what you are really interested in. Occasionally, though, you'll need to actually calculate the damping coefficient and determine whether the system is overdamped or not.

The first time that I ran into this problem was in an automatic controls class where I needed to calculate the damping coefficient for an elastic band in order to use it elsewhere in our equations. At that time I only needed to calculate the coefficients for a few different configurations so I did the analysis in Excel. Essentially what you do is use the solver plug-in to minimize the sum of the error squared between your theoretical data and the experimental data at the peaks of the wave. It works pretty well and here's an [ example] file to see what I mean.

The next time I ran into this problem I was working on a camshaft lab where we were analyzing the valve spring behavior for a range of springs, shims and engine speeds. We were generating a crap load of data, something like a few million data points in all. Using Excel wasn't very practical in this case, so I wrote a nice automated Matlab script to do it for me.

We're assuming that the amplitude envelope is described by a decaying exponential function. The problem is that y=a*exp(sigma*x) isn't particularly linear. To get around this problem we take the log of both sides to get log(y) = log(a) + sigma*x. On a quick side note, log means the natural logarithm (ln) and b will be negative because we are dealing with decay.

Now that we have a linear equation we can easily do a least-squares error fit in order to calculate the constants a and sigma. We'll optimize the equation Y = A + BX where Y = log(y), A = log(a), B = sigma and X = x. After some linear algebra you can show that p = (A'*A)^-1*A*Y where A = [1 X1; 1 X2 ... 1 Xn], Y = [Y1; Y2 ... Yn] and p = [a;sigma]. In Matlab you can just use p = A\Y. You can calculate a and b from A and B and then plot away. But what is the physical meaning of sigma and what does it tell you about the damping for the oscillator?

The idea is that sigma = chi*omega_n where chi is the damping coefficient and omega_n is the natural frequency. The damped natural frequency is equal to sqrt(omega_n^2-sigma^2). There are three cases then:

  • Underdamped: omega_n^2 - sigma^2 > 0
  • Critically damped: omega_n^2 - sigma^2 = 0
  • Overdamped: omega_n^2 - sigma^2 < 0

Here's some example code that generates sample data and calculates the damping parameters.

clf
clear all

% This is just example data
xData = linspace(0,20*pi,1000);
yData = cos(2*xData).*exp(-.05*xData);

% If we are only interested in a certain data range
rangeLow  = 0;
rangeHigh = 10000;

xCoords = [];
yCoords = [];

% Find the x,y coordinates for the oscillation peaks
yData    = yData(:);
upOrDown = sign(diff(yData));
maxFlags = [upOrDown(1)<0 ; diff(upOrDown)<0 ; upOrDown(end)>0];
maxIndices   = find(maxFlags);

for ii = 1:length(maxIndices)
    if (maxIndices(ii) > rangeLow) && (maxIndices(ii) < rangeHigh)
        xCoords = [xCoords xData(maxIndices(ii))];
        yCoords = [yCoords yData(maxIndices(ii))];
    end
end

% Calculate the optimal values of a and b
A = ones(length(xCoords),2);
Y = ones(length(xCoords),1);

% We take the log of the actual data (yCoords)
for ii = 1:length(xCoords)
    A(ii,2) = xCoords(ii);
    Y(ii,1) = log(yCoords(ii));
end

% After finding A and B, we know that a = exp(A) and b = B
p = A\Y;
a = exp(p(1));
b = p(2);

% We plug them into their function and plot away. Note that b < 0
curvePlot = a*exp(b*xData);

figure(1)
plot(xData,yData,'-b',xCoords,yCoords,'or',xData,curvePlot,'--g');
legend('Data','Points','Curve');

The problem with this code, however, is that if your data has noise in it then it will have trouble finding the actual oscillator peaks. Currently the code considers a peak to a point where the first derivative is positive before and negative after a point. We can quickly improve on this code by incorporating a tolerance band and iterating through the data to find maxima. The tolerance band specifies how far on either side of a point to check in order to verify that the current point is the maximum value in that range. Iterating is a bit slower than the 'find' approach, but it'll do.

The following code does a lot more than the previous code. It was used in analyzing the damping of valve spring oscillations at high engine speeds, so it:

  • Loads the spring force data for a specified engine speed. Three springs and three shim sizes were studied.
  • Shifts the force data so that the average force during oscillations is zero. This is necessary accurately report the magnitude of oscillation, which is handy for diagnosing dynamic coil binding.
  • Calculates the damping factor.
  • Calculates the theoretical natural frequency from known masses and spring constants.
  • Calculates the theoretical damped natural frequency, which can later be compared to a fast Fourier transform (FFT) of the system to verify the numbers.
clf
clear all

engineSpeed = 4800;

temp = load([num2str(engineSpeed) '_Ford_Damped_NoShim.m']);
AllData(:,1) = temp(:,2);
temp = load([num2str(engineSpeed) '_Ford_Damped_Shim034.m']);
AllData(:,2) = temp(:,2);
temp = load([num2str(engineSpeed) '_Ford_Damped_Shim064.m']);
AllData(:,3) = temp(:,2);

temp = load([num2str(engineSpeed) '_CASS_Damped_NoShim.m']);
AllData(:,4) = temp(:,2);
temp = load([num2str(engineSpeed) '_CASS_Damped_Shim034.m']);
AllData(:,5) = temp(:,2);
temp = load([num2str(engineSpeed) '_CASS_Damped_Shim064.m']);
AllData(:,6) = temp(:,2);

temp = load([num2str(engineSpeed) '_CADS_NoShim.m']);
AllData(:,7) = temp(:,2);
temp = load([num2str(engineSpeed) '_CADS_Shim034.m']);
AllData(:,8) = temp(:,2);
temp = load([num2str(engineSpeed) '_CADS_Shim064.m']);
AllData(:,9) = temp(:,2);

for ii =  1:9
    xData = linspace(0,720,2048)/engineSpeed/360*60;
    yData = AllData(:,ii);
    yData = yData - mean(yData(100:500));

    % If we are only interested in a certain data range
    rangeLow  = 101;
    rangeHigh = 1200;
    checkRange = 100;

    maxIndices = [];

    for jj = 1:length(yData)
        if (jj > rangeLow) & (jj < rangeHigh)
            for kk = -checkRange:checkRange
                if (yData(jj) >= yData(jj+kk))
                    if kk == checkRange
                        maxIndices = [maxIndices jj];
                    end
                else
                    break
                end
            end
        end
    end

    xCoords = xData(maxIndices);
    yCoords = yData(maxIndices);

    % Calculate the optimal values of a and b
    A = ones(length(xCoords),2);
    Y = ones(length(xCoords),1);

    % We take the log of the actual data (yCoords)
    for jj = 1:length(xCoords)
        A(jj,2) = xCoords(jj);
        Y(jj,1) = log(yCoords(jj));
    end

    % After finding A and B, we know that a = exp(A) and b = B
    p = A\Y;
    a(ii) = exp(p(1));
    sigma(ii) = p(2);

    % We plug them into their function and plot away. Note that b < 0
    curvePlot = a(ii)*exp(sigma(ii)*xData);
    figure(ii)
    plot(xData,yData,'-b',xCoords,yCoords,'or',xData,curvePlot,'--g');
    title('CADS (4800 RPM)','FontSize',20);
    xlabel('Time (s)','FontSize',20);
    ylabel('Spring Force (N)','FontSize',20);
    legend('Experimental Data','Analyzed Points','Fitted Curve');
    pause
    print -r300 -djpeg80 SpringOscillations
end

% Calculate the natural frequencies:
mretainer = 28e-3; %kg
mkeepers  = 4.3e-3;
mrocker   = 116.5e-3;
mvalve    = 96.4e-3;

msys = mretainer + mkeepers + mrocker + mvalve;

mspring(1) = 61.6e-3; % Ford
mspring(2) = mspring(1);
mspring(3) = mspring(1);
mspring(4) = 75.4e-3; % CASS
mspring(5) = mspring(4);
mspring(6) = mspring(4);
mspring(7) = 70.2e-3; % CADS
mspring(8) = mspring(7);
mspring(9) = mspring(7);

meff = mspring*.33 + msys;

k(1) = 28.3e3; % Damped Ford (N/m)
k(2) = k(1);
k(3) = k(1);
k(4) = 50.4e3; % Damped CASS
k(5) = k(4);
k(6) = k(4);
k(7) = 31.4e3; % Damped CADS
k(8) = k(7);
k(9) = k(7);

w_natural = sqrt(k./meff)
w_damped = sqrt(w_natural.^2-sigma.^2)

Note that I wrote all of this code way back in 2005 or 2006, and since then have become a much better Matlab coder. So take the code in this post with a grain of salt.

Run Comsol Simulations from Matlab

Thursday, March 1st, 2007

Note: This post is slightly outdated. Since I wrote this post back in 2004, Comsol integrated parameter looping to their solver (it was called Femlab back then). Many use cases that would have previously required Matlab integration can be handled entirely within Comsol now. The particular case that is shown could not (because the permittivity is calculated from a helper function). Unfortunately I do not have the time to answer specific questions or basic FEA issues, but hope that this post gets you pointed in the right direction for simulation automation.

Comsol is extremely powerful because it's well integrated with Matlab. However, the ability to build and run simulations in Matlab isn't well documented.

Here are a few cases where you would want to build a simulation in Matlab:

  • You want to simulate a variety of geometries to optimize a system
  • You are generating tons of data and you would like to save the results for later analysis and plot generation
  • Your simulation is so complicated that the file is becoming tens or hundreds of megabytes in size

Instructions

Create a new simulation in Comsol or open an existing one. Now, choose Save As... from the File menu. Before you do anything else, look at the options that you have for saving it. One of them is to save it as a .m file. Do that and then open up the existing file in a text editor. You'll see something like this:

flclear fem
clear vrsn
vrsn.name = 'FEMLAB 3.1';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 157;
vrsn.rcs = '$Name: $';
vrsn.date = '$Date: 2004/11/12 07:39:54 $';
fem.version = vrsn;
...

Take a few minutes to dig through the file and guess what it's doing. It's important to realize that this is a normal .m file so you have access to all of Matlab's libraries and commands.

So how do you run these things? The next time that you start up Comsol be sure to use the Comsol with Matlab option and then just open the file in Matlab like you would any other script. This last step gave me plenty of grief and it's surprising that Comsol doesn't more clearly communicate the massive power that you gain by working directly with the text file (edit: at least in 2004 when I was working on this, I've heard rumors that things are a lot better documented now).

My usual procedure for using this feature is create a new simulation in Comsol using the GUI with the exact options, geometry, boundary conditions, etc. that I want. Once it is done, run the simulation and then save it as a .m file. This approach is useful because Comsol will generate all of the hard part for you and you can just go back and add your for loops, etc.

Here's an example file that I used for simulating the local electromagnetic field enhancement for a gold nanocrescent at UC Berkeley.

% Joseph Doll
% OpticalOneLayer.m
% FEMLAB 3.1 simulation to calculate the maximum local electric field
% enhancement for a gold or silver nanocrescent of arbitrary ID, OD,
% and opening width immersed in a low absorption dielectric medium of
% arbitrary dielectric constant.
function [simResults, lambda] = NanocrescentSim(OD, IDRatio, ...
  DRatio, lambda, eps_o)

flclear fem
clear vrsn
vrsn.name = 'FEMLAB 3.1';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 157;
vrsn.rcs = '$Name: $';
vrsn.date = '$Date: 2004/11/12 07:39:54 $';
fem.version = vrsn;

% -------------------------------------------------------------
% -------------- Start Geometry Initialization ----------------
% -------------------------------------------------------------

% Fix outside diameter
% Specify inside diameter as a proportion of OD
% Specify delta as a proportion of ID

indexCounter = 1;
theta = 0;

for outsideRadius = OD/2
  for insideRadius = outsideRadius.*IDRatio
    for delta = (insideRadius*2).*DRatio

    sprintf('OD = %d, ID = %d, OW = %d', ...
      2*outsideRadius, 2*insideRadius, delta)

    % Calculate the horizontal displacement of each ellipse based
    % upon the dimensions noted above (radii and opening widths)
    outerCircleDx = 0;
    outerCircleDy = 0;

    innerCircleDx = -((outsideRadius - insideRadius) + ...
      (insideRadius*(1-cos(asin(delta/(2*insideRadius))))) - ...
      (outsideRadius*(1-cos(asin(delta/(2*outsideRadius))))));
    innerCircleDy = 0;

    % Define the ellipses based upon the calculated dimensions
    g1=ellip2(outsideRadius,outsideRadius,'base','center','pos', ...
      [outerCircleDx,outerCircleDy],'rot','0');
    g5=ellip2(insideRadius,insideRadius,'base','center','pos', ...
      [innerCircleDx,innerCircleDy],'rot','0');

    % Define the environment bound rectangle
    bound_rect = rect2('4e-6','2e-6','base','center', ...
      'pos',{'0','0'},'rot','0');

    % Boolean operations to create the "sharp structures"
    sharp_nc = geomcomp({g1,g5},'ns',{'g1','g5'},'sf', ...
      'g1-g5','edge','none');

    % Fillets
    round_nc = fillet(sharp_nc,'radii',2E-9,'point',[1 2]);

    % Rotate the finished nanocresent
    nanocrescent = rotate(round_nc, theta ,[0,0]);

    clear s
    s.objs = {nanocrescent, bound_rect};
    s.name={'NC','BR'};
    s.tags={'nanocrescent','boundrect'};

    fem.draw=struct('s',s);
    fem.geom=geomcsg(fem);

    % ---------------------------------------------------------
    % ---------------- Start Simulation Loop ------------------
    % ---------------------------------------------------------

    % Vector of the max electric field at each wavelength
    emax = [];

    % Calculate the dielectric constant values
    % Requires helper function OpticalData.m
    eps_i = OpticalData( lambda(1), ...
      lambda(length(lambda)),length(lambda),1);

    for ii = 1:length(lambda)
      fem.const={'eps_i',eps_i(ii),'eps_o',eps_o};

      clear appl
      appl.mode.class = 'InPlaneWaves';
      appl.module = 'CEM';
      appl.assignsuffix = '_wh';

      % Solve via input wavelength
      clear prop
      prop.field='TM';
      prop.inputvar='lambda';
      appl.prop = prop;

      % Set the boundary conditions
      clear bnd
      bnd.H0 = {{0;0;1},{0;0;0},{0;0;0}};
      bnd.type = {'NR','cont','NR'};
      bnd.ind = [1,3,3,3,2,2,2,2,2,2,2,2,2,2];
      appl.bnd = bnd;

      % Set the domain variables
      clear equ
      equ.epsilonr = {'eps_o','eps_i'};
      equ.ind = [1,2];
      appl.equ = equ;

      % Set the current wavelength
      appl.var = {'nu','1e9','lambda0',lambda(ii)};
      fem.appl{1} = appl;
      fem.border = 1;

      fem=multiphysics(fem);

      % Mesh initialize parameters, smaller = finer
      fem.mesh=meshinit(fem, ...
        'hmaxfact',1, ...
        'hcurve',0.1, ...
        'hgrad',1.2, ...
        'hcutoff',0.001);

      fem.xmesh=meshextend(fem);

      fem.sol=femlin(fem, ...
        'solcomp',{'Hz'}, ...
        'outcomp',{'Hz'});

      % posteval returns the electric field at each point
      eii_struct = posteval(fem,'normE_wh','dl',[1]);

      % Take the maximum field value and add it to emax
      eii_max = max(eii_struct.d);
      emax = [emax eii_max];
    end

    simResults{indexCounter} = emax;
    indexCounter = indexCounter + 1;
  end
end
end

Using this approach I was able to automate process of finding the peak electric field intensity as a function of design parameters. Combined with another script I was able to run a simulation for about 3 days on end that tested 300+ different geometry combinations. And all I had to do was set it up, click a button, and put a note on the computer to not disturb it. Like I said, pretty awesome.


fetishroom