光学薄膜設計ソフトの設計 その28 - 材料屈折率データベース [考え中 - 光学薄膜設計]
何回か前、材料の屈折率の最適化に関する問題を考えた。最適化の前にまず、実在する材料の屈折率/分散のデータが必要になる。
TFCalcのサイトで見つけたSopra SA社が公開しているデータを取り込むことを考えた。そのやり方を整理する。
Sopra SAのデータから、
- 誘電体はSellmeierの分散式
- 金属はDrudeの分散式（複数のプラズマ振動の和）
誘電体は線形化したSellmeierの式
しょうがないので、虚数部を持っているデータは、
- 一旦全部を3次のスプラインで内挿
- 波長で均一な刻みの表を作り直す
- 表から線形内挿で波長に対する値を決める
最終的にはメインバンドルのリソースとしてアプリケーションに含めて、計算エンジンから要求されたらデータを読み込む、というようにしたい。そのためにはプラグインの形式にする、というのでもいいけど、もっと簡単に、データだけ書いたファイルを作って、必要になったとき読み込むようにする。データはテキストファイルにしておく。
どうせ自分で使うだけなので、えいや、でデータ用のテキストファイルの形式を決めてしまった。
| 1行目 | タイプ | 整数 |
| 2行目 | 説明 | 文字列 |
| 3行目 | 最短波長 | 実数 |
| 4行目 | 最長波長 | 実数 |
- Sellmeier
- Drude
- 表
- その他
| 5行目 | 波長刻み | 実数 |
| 6行目 | 波長、屈折率の実数部、屈折率の虚数部 | 3つの実数 |
とする。
これを読み込んで波長の値に対して屈折率の値を返すクラスを作る。OTFMediumModelとしよう。
#import "OTFComplex.h"
extern NSString *nullModelName;
typedef struct refractiveIndexCoefficient *riPtr;
@interface OTFMediumModel : NSObject {
NSString *modelName;
NSString *descr;
riPtr riCoef;
double complex (*refFunc)();
}
+ (NSArray *)modelNames;
- (id)initWithModelName:(NSString *)mName;
- (NSString *)modelName;
- (double complex)refractiveIndexAt:(double)lambda;
// Private methods
- (void)loadCoefficient;
- (void)setSellmeierFrom:(NSArray *)lines;
- (void)setTableFrom:(NSArray *)lines;
@end
このクラスはinitializeクラスメソッドでさっきのファイルの入っているディレクトリを見て、ファイルの名前をmodelNamesに設定する。
static NSMutableArray *modelNames = nil;
static NSMutableArray *modelPaths;
static NSString *modelExtension = @"ind";
+ (void)initialize
{
if (modelNames)
return;
NSArray *paths = [[NSFileManager defaultManager]
directoryContentsAtPath:modelDirectoryPath];
int tmpCount = [paths count];
NSEnumerator *en = [paths objectEnumerator];
id obj;
modelNames = [[NSMutableArray alloc] initWithCapacity:tmpCount];
modelPaths = [[NSMutableArray alloc] initWithCapacity:tmpCount];
while (obj = [en nextObject]) {
if ([[obj pathExtension] isEqualToString:modelExtension]) {
id path = [modelDirectoryPath stringByAppendingPathComponent:obj];
id name = [[obj lastPathComponent] stringByDeletingPathExtension];
[modelPaths addObject:path];
[modelNames addObject:name];
}
}
if (![modelNames containsObject:nullModelName])
[NSException raise:@"Defective Database" format:@"Lack of Null Medium"];
allocatedObjects = [[NSMutableSet alloc] initWithCapacity:0];
}
+ (NSArray *)modelNames
{
return modelNames;
}
クラス変数のmodelNamesとmodelPathsに拡張子がindのものを選んで配列に入れているだけ。クラスメソッドのmodelNamesはその配列を返している。allocatedObjectsは
static NSMutableSet *allocatedObjects;で、OTFMergableと同じように、同一のモデル名を持つインスタンスはひとつだけに限定するためのもの。OTFMergableとは違う理由なのでOTFMergableが親クラスにはなっていない。もしOTFMergableを抽象クラスではなく、プロトコルにするならそれをadoptするほうがいいだろう。
インスタンス生成のためのinit...メソッドは
- (id)initWithModelName:(NSString *)mName
{
self = [super init];
id obj;
unsigned int idx = [modelNames indexOfObject:mName];
if (idx == NSNotFound) {
[self release];
return nil;
}
modelName = [modelNames objectAtIndex:idx]; // (1)
riCoef = NULL; // flag for local allocation
if (obj = [allocatedObjects member:self]) {
[self release];
return [obj retain];
}
riCoef = (refractiveIndexCoefficient *)
malloc(sizeof(refractiveIndexCoefficient));
[self loadCoefficient];
[allocatedObjects addObject:self];
return self;
}
渡された名前がmodelNames配列の中になければデータがないと言うことなのでnilを返す。また、すでにallocatedObjectsの中に同じ名前のオブジェクトがあれば自分を解放してそちらを返す。
NSSetの中で同じかどうかは昔試したように、hashとisEqual:メソッドで比較されるので
- (BOOL)isEqual:(id)obj
{
return modelName == [obj modelName];
}
- (NSUInteger)hash
{
return (NSUInteger)modelName;
}
と書いておいて名前が同じものをさしている場合だけ同じと見なす。これは(1)のところで、渡されてきた文字列ではなくてmodelNames配列から取り出しているので文字列比較でなくてもかまわない。この場合はhashで比較するのとisEqual:での結果が一致する。でもこういうのはあまり美しくないな、速いけど。
同じのがなければ係数の領域を確保してファイルの中身を使って初期化する。係数の初期化はloadCoefficientメソッドの中で行う。
- (void)loadCoefficient
{
NSError *error;
NSString *path = [modelPaths objectAtIndex:
[modelNames indexOfObject:modelName]];
NSString *coefString
= [NSString stringWithContentsOfFile:path
encoding:NSASCIIStringEncoding
error:&error];
NSArray *lines = [coefString arrayOfLines];
switch ([[lines objectAtIndex:0] intValue]) {
case OTFDielectricMediumType:
[self setSellmeierFrom:lines];
break;
case OTFTableMediumType:
[self setTableFrom:lines];
break;
default:
break;
}
}
ファイルパスを設定して中身をNSStringに読み込む。それを一行ごとに分ける。その1行目を読んでSellmeierとしてセットするのか、表とするのかを決める。1行ごとに分けるのはHMDTの「1行ずつSubstringを取り出す」をそのままいただきました。
屈折率係数の構造体は
typedef struct refractiveIndexCoefficient {
int typ;
spectralRange rng;
union {
sellmeierCoefficients sel;
drudeCoefficients dru;
tableCoefficients tbl;
} coef;
} refractiveIndexCoefficient;
で、typは係数のタイプを表す。つぎのrngは
typedef struct {
double min;
double max;
} spectralRange;
で、波長の適応範囲を保持する。
つぎは、タイプによって違う形式の構造体のunionになっている。
Sellmeierの場合、
typedef struct sellmeierCoefficients {
double a;
double b;
double c;
double d;
double e;
double f;
double g;
} sellmeierCoefficients;
で、係数を並べただけ。
表形式の場合も
typedef struct {
double step;
unsigned int len;
double complex *ri;
} tableCoefficients;
刻みの幅と長さと値の配列を持っているだけ。そのまんま。
例えばSellmeierはどう読み込んでいるかと言うと
- (void)setSellmeierFrom:(NSArray *)lines
{
NSEnumerator *en = [lines objectEnumerator];
refFunc = refractiveIndexFromSellmeier;
[en nextObject]; // type
descr = [[en nextObject] retain]; // description
riCoef->typ = OTFDielectricMediumType;
riCoef->rng.min = [[en nextObject] doubleValue];
riCoef->rng.max = [[en nextObject] doubleValue];
riCoef->coef.sel.a = [[en nextObject] doubleValue];
riCoef->coef.sel.b = [[en nextObject] doubleValue];
riCoef->coef.sel.c = [[en nextObject] doubleValue];
riCoef->coef.sel.d = [[en nextObject] doubleValue];
riCoef->coef.sel.e = [[en nextObject] doubleValue];
riCoef->coef.sel.f = [[en nextObject] doubleValue];
riCoef->coef.sel.g = [[en nextObject] doubleValue];
}
もうそのまんま。表はもういいでしょう、同じことだから。
ここで設定されているrefFuncというインスタンス変数はC関数へのポインタでSellmeierの場合、refractiveIndexFromSellmeierという関数をさすようにしている。この関数は
static double complex refractiveIndexFromSellmeier(
refractiveIndexCoefficient *coefs,
double lambda)
{
if (lambda < coefs->rng.min)
lambda = coefs->rng.min;
else if (lambda > coefs->rng.max)
lambda = coefs->rng.max;
double l2 = lambda * lambda;
double invl2 = 1.0 / l2;
sellmeierCoefficients *sel = &(coefs->coef.sel);
return (double complex)(sqrt(sel->a
+ l2 * (sel->b + sel->g * l2)
+ invl2 * (sel->c
+ invl2 * (sel->d
+ invl2 * (sel->e
+ invl2 * sel->f)))));
}
のように係数へのポインタと波長の値を引数にもらって屈折率の値を返す関数。
波長が定義されている範囲を超えているとき、そのまま計算すると変な値が返るので、範囲の端の値を返すことにする。
オブジェクトに対して屈折率の値の要求があったとき、これで計算する。
- (double complex)refractiveIndexAt:(double)lambda
{
return (*refFunc)(riCoef, lambda);
}
まあ、どうということもない。
表で屈折率が表されている方も、面倒なだけで同じことなので省略する。
屈折率モデルにする形式が増えたら、この格好で増やしていけばいい。
ちなみに、例えばSiO2（二酸化硅素）のファイル「SIO2.ind」は
1 SiO2-Silicon Dioxide 0.163138396 2.06641969 2.1033415981683974 -0.008488713592969977 0.008926516328328488 0.00010478955297666573 -2.340821478762694e-6 1.3595019653044792e-7 -0.0002254782075299302で、
OTFMediumModel *mmodel = [[OTFMediumModel alloc] initWithModelName:@"SIO2"];
for (l = 0.2 ; l <= 0.8 ; l += 0.01)
printf("\%f\t\%f\n", l, creal([mmodel refractiveIndexAt:l]));
として出力したものをMathematicaでプロットしてみると図-5.6となる。
現実のものとあっているのかどうかわからんけど、いかにもそれっぽい。
今日はなんかごちゃごちゃ細かく書きすぎた。こんなのSourceをみればわかるな。
2009-02-03 23:58
nice!(0)
コメント(0)
トラックバック(0)
コメント 0