From 7665696e387c447d59ecb80b2e14705ab1c61568 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Wed, 22 Jan 2014 14:50:13 +0000 Subject: [PATCH] added rge-averaging support --- .../libPhotoMeissner/test/photo_meissner.cpp | 860 ++++++++++++++---- 1 file changed, 705 insertions(+), 155 deletions(-) diff --git a/src/external/libPhotoMeissner/test/photo_meissner.cpp b/src/external/libPhotoMeissner/test/photo_meissner.cpp index 1c5bf8fd..f9838d74 100644 --- a/src/external/libPhotoMeissner/test/photo_meissner.cpp +++ b/src/external/libPhotoMeissner/test/photo_meissner.cpp @@ -40,12 +40,32 @@ using namespace std; #include +typedef struct { + double lamL; + double lamP; + double z0; + double deadLayer; + double filmThickness; +} PhotoParam; + +double InuMinus(double nu, double x); + //-------------------------------------------------------------------------- void photo_meissner_syntax() { - cout << endl << "usage: photo_meissner lamL lamP z0 deadLayer [filmThickness]"; - cout << endl << " {--at | --range ]} [--dump ] | "; - cout << endl << " --version | [--help]" << endl; + cout << endl << "usage: (0) photo_meissner --version | [--help]" << endl; + cout << endl << " (1) photo_meissner lamL lamP z0 deadLayer [filmThickness]"; + cout << endl << " {--at | --range ]} [--dump ]" << endl; + cout << endl << " this option is used to calculate fields without any muon stopping"; + cout << endl << " distribution averaging" << endl; + cout << endl << " (2) photo_meissner lamL lamP z0 deadLayer [filmThickness] --file "; + cout << endl << " [--dump ]" << endl; + cout << endl << " this option is used to calculate and , for a given muon"; + cout << endl << " implantation energy, given by the rge-file." << endl; + cout << endl << " (3) photo_meissner lamL lamP z0 deadLayer [filmThickness]"; + cout << endl << " --prefix --energies [--dump ]" << endl; + cout << endl << " this option is used to calculate and , for a whole list of muon energies"; + cout << endl << " as given in the rge-file list which will be generated by and ." << endl; cout << endl << " lamL : London penetration depth in (nm)"; cout << endl << " lamP : Photo induced penetration depth (nm)"; cout << endl << " z0 : length scale of the photo induced effect (nm)"; @@ -55,12 +75,489 @@ void photo_meissner_syntax() cout << endl << " --at : calculates the field at the position (nm)."; cout << endl << " --range : calculate the field in the range from"; cout << endl << " to in (nm)."; - cout << endl << " either --at or --range needs to be present."; + cout << endl << " either --at or --range needs to be present."; cout << endl << " --dump : dump the result into the file ."; + cout << endl << " if this option is not given, the result is dumped to stdout"; + cout << endl << " --file : muon stopping distribution to calculate ."; + cout << endl << " --prefix : prefix to generate the rge-path-file-name, e.g. --prefix YBCO_E"; + cout << endl << " --energies : energy list, where list is a comma seprated energy list in eV,"; + cout << endl << " e.g. --energies 1000,1500,2000,2500,...,250000. Optionally, it can also"; + cout << endl << " be --energies 1000:25000:1000, in which case it means start at 1000eV"; + cout << endl << " increment by 1000eV, and fill the energy list up to 25000eV."; + cout << endl << " NO spaces between the comma separated list allowed!"; cout << endl << endl; } //-------------------------------------------------------------------------- +/** + *

Routine to tokenize a string. Empty tokens will be ignored. + * + * \param cstr string to be tokenized + * \param delim string of characters to be parsed + * \param tok vector of tokens + */ +void tokenizer(const char *cstr, const char *delim, vector &tok) +{ + tok.clear(); + char sstr[128]; + + int ss=0; + for (int i=0; i0) // ignore empty tokens + tok.push_back(sstr); + ss = i+1; + } + } + } + + // handle last token if the last character is not a delim + for (int i=0; i0) + tok.push_back(sstr); + } + } +} + +//-------------------------------------------------------------------------- +/** + *

+ * + * \param prefix used to generate the path file name + * \param elist is the list of energies, which is either a comma sparated list, or has the format :: + * \param rgeFln rge-file list + */ +int generateRgeFln(const string prefix, const string elist, vector &rgeFln) +{ + vector tok; + string str(""); + char istr[128]; + int start, end, step, ival; + + if (elist.find(",") != string::npos) { // comma separated energy list + tokenizer(elist.c_str(), ",", tok); + for (unsigned int i=0; i::. Found " << tok.size() << " tokens instead of 3." << endl; + return 0; + } + // check that all tokens are indeed numbers > 0 + start = atoi(tok[0].c_str()); + end = atoi(tok[1].c_str()); + step = atoi(tok[2].c_str()); + if ((start <= 0) || (end <= 0) || (step <= 0) || (end <= start)) { + cerr << endl << "**ERROR** in generateRgeFln: option --energies ::. One of these value is not a number > 0." << endl; + return 0; + } + ival = start; + do { + sprintf(istr, "%d", ival); + str = prefix + istr + ".rge"; + rgeFln.push_back(str); + ival += step; + } while (ival <= end); + } else { // something is wrong + cerr << endl << "**ERROR** in generateRgeFln: found wrong format for option --energies: " << elist << endl; + return 0; + } + return 1; +} + +//-------------------------------------------------------------------------- +/** + *

+ * + * \param line of a rge-file. + * \param z value: position in nm + * \param n value: number of particles at position z + */ +int getRgeValue(const char *line, double &z, double &n) +{ + vector tok; + tokenizer(line, " ", tok); + + z = -1.0; n = -1.0; + if (tok.size() == 2) { + if (isdigit(tok[0][0])) { + z = strtod(tok[0].c_str(), 0); + } else { + tok.clear(); + return 1; + } + if (isdigit(tok[1][0])) { + n = strtod(tok[1].c_str(), 0); + } else { + tok.clear(); + return 1; + } + } else { + tok.clear(); + return 1; + } + tok.clear(); + + if ((z == -1.0) || (n == -1.0)) { + cerr << endl << "**ERROR in getRgeValue(): line '" << line << "' found. Seems to be invalid." << endl; + return 0; + } + + return 1; +} + +//-------------------------------------------------------------------------- +/** + *

+ * + * \param fln rge-file name + * \param z position vector in nm + * \param n number of particle vector corresponding to z + */ +int readRgeFile(const char *fln, vector &z, vector &n) { + + ifstream fin(fln, ifstream::in); + if (!fin.is_open()) { + cerr << endl << "**ERROR** couldn't open '" << fln << "', will quit." << endl; + return 0; + } + + char line[512]; + double zz, nn; + z.clear(); + n.clear(); + while (fin.good()) { + fin.getline(line, 512); + if (!getRgeValue(line, zz, nn)) { + fin.close(); + return 0; + } + if (zz != -1.0) { // no comment or empty line + z.push_back(zz/10.0); // A -> nm + n.push_back(nn); + } + } + + fin.close(); + + return 1; +} + +//-------------------------------------------------------------------------- +/** + *

+ * + * \param fln rge-path-file-name vector + * \param zz all position vectors in nm + * \param nn all number of particle vectors corresponding to zz + */ +int readAllRgeFiles(const vector &fln, vector< vector > &zz, vector< vector > &nn) +{ + for (unsigned int i=0; i get interpolated number of particles at position z for a given rge-set. + * + * \param z requested position (in nm). + * \param zz position vector (in nm). + * \param nn number of particle vector corresponding to zz + */ +double get_n(const double z, vector &zz, vector &nn) +{ + int idx=-1; + double n=0.0; + + // find proper position index + for (unsigned int i=0; i z) { + idx = (int)i-1; + break; + } + } + + if (idx == -1) + return 0.0; + + // linear interpolation + n = nn[idx] + (nn[idx+1]-nn[idx])*(zz[idx+1]-z)/(zz[idx+1]-zz[idx]); + + return n; +} + +//-------------------------------------------------------------------------- +/** + *

interpolated number of particle distribution. + * + * \param rgeZ position values of an rge-vector (in nm) + * \param rgeN number of particle vector corresponding to rgeZ + * \param zz interpolated position vector + * \param nn interpolated number of particle vector corresponding to zz + */ +void matched_nz(vector &rgeZ, vector &rgeN, vector &zz, vector &nn) +{ + for (unsigned int i=0; iCalculates photo Meissner screening field at position z for a given set of parameters param. + * + * \param z position at which the screening profile shall be calculated (in nm). + * \param param structure containing all relvant parameters + */ +double calcSingleField(const double z, const PhotoParam param) +{ + // estimate dz + double dz = 0.0; + if (param.lamL < param.lamP) + dz = param.lamL / 500.0; + else + dz = param.lamP / 500.0; + + double nuPlus=0.0; + double nuPhoto=0.0; + double b=-1.0, b1=0.0; + if (param.filmThickness == -1.0) { // semi-infinite + nuPlus = 2.0*param.z0/param.lamL; + nuPhoto = 2.0*param.z0/param.lamP; + if (z < param.deadLayer) { + b = 1.0; + } else { + b1 = gsl_sf_bessel_Inu(nuPlus, nuPhoto); + assert(b1>0.0); + b = gsl_sf_bessel_Inu(nuPlus, nuPhoto * sqrt(exp(-(z-param.deadLayer)/param.z0)))/b1; + } + } else { // film + nuPlus = 2.0*param.z0/param.lamL; + nuPhoto = 2.0*param.z0/param.lamP; + double beta = exp(-(param.filmThickness/2.0-param.deadLayer)/param.z0); + double NN = InuMinus(nuPlus, nuPhoto*beta)*gsl_sf_bessel_Inu(nuPlus, nuPhoto) - + InuMinus(nuPlus, nuPhoto)*gsl_sf_bessel_Inu(nuPlus, nuPhoto*beta); + assert(NN>0); + if ((z < param.deadLayer) || (z > param.filmThickness-param.deadLayer)) { + b = 1.0; + } else { + b = (InuMinus(nuPlus, nuPhoto*exp(-(z-param.deadLayer)/(2.0*param.z0)))*(gsl_sf_bessel_Inu(nuPlus, nuPhoto)-gsl_sf_bessel_Inu(nuPlus, nuPhoto*beta))- + gsl_sf_bessel_Inu(nuPlus, nuPhoto*exp(-(z-param.deadLayer)/(2.0*param.z0)))*(InuMinus(nuPlus, nuPhoto)-InuMinus(nuPlus, nuPhoto*beta)))/NN; + } + } + return b; +} + +//-------------------------------------------------------------------------- +/** + *

Calculates photo Meissner screening field between zStart and zEnd for a given set of parameters param. + * + * \param zStart start position from where to calculate B(z) + * \param zEnd end position up to where to calculate B(z) + * \param zz position vector in nm (return value) + * \param bb photo Meissner screening vector corresponding to zz + * \param bbL London Meissner screening vector corresponding to zz + */ +void calcFields(const double zStart, const double zEnd, const PhotoParam param, vector &zz, vector &bb, vector &bbL) +{ + // estimate dz + double dz = 0.0; + if (param.lamL < param.lamP) + dz = param.lamL / 500.0; + else + dz = param.lamP / 500.0; + + double nuPlus = 2.0*param.z0/param.lamL; + double nuPhoto = 2.0*param.z0/param.lamP; + + double z = zStart; + double b=-1.0, bL=-1.0; + if (param.filmThickness == -1.0) { // semi-infinite + double b1 = gsl_sf_bessel_Inu(nuPlus, nuPhoto); + assert(b1>0.0); + do { + if (z < param.deadLayer) { + b = 1.0; + bL = 1.0; + } else { + b = gsl_sf_bessel_Inu(nuPlus, nuPhoto * sqrt(exp(-(z-param.deadLayer)/param.z0)))/b1; + bL = exp(-(z-param.deadLayer)/param.lamL); + } + zz.push_back(z); + bb.push_back(b); + bbL.push_back(bL); + z += dz; + } while (z <= zEnd); + } else { // film + double beta = exp(-(param.filmThickness/2.0-param.deadLayer)/param.z0); + double NN = InuMinus(nuPlus, nuPhoto*beta)*gsl_sf_bessel_Inu(nuPlus, nuPhoto) - + InuMinus(nuPlus, nuPhoto)*gsl_sf_bessel_Inu(nuPlus, nuPhoto*beta); + assert(NN>0); + do { + if ((z < param.deadLayer) || (z > param.filmThickness-param.deadLayer)) { + b = 1.0; + bL = 1.0; + } else { + b = (InuMinus(nuPlus, nuPhoto*exp(-(z-param.deadLayer)/(2.0*param.z0)))*(gsl_sf_bessel_Inu(nuPlus, nuPhoto)-gsl_sf_bessel_Inu(nuPlus, nuPhoto*beta))- + gsl_sf_bessel_Inu(nuPlus, nuPhoto*exp(-(z-param.deadLayer)/(2.0*param.z0)))*(InuMinus(nuPlus, nuPhoto)-InuMinus(nuPlus, nuPhoto*beta)))/NN; + bL = cosh((param.filmThickness/2.0-z)/param.lamL)/cosh((param.filmThickness/2.0-param.deadLayer)/param.lamL); + } + zz.push_back(z); + bb.push_back(b); + bbL.push_back(bL); + z += dz; + } while (z <= zEnd); + } +} + +//-------------------------------------------------------------------------- +/** + *

calculates the average of a vector aa based on the distribution nn. + * + * \param aa vector to be averaged + * \param nn distribution over which aa has to be averaged + */ +double averaged(vector &aa, vector &nn) +{ + if (aa.size() != nn.size()) { + cerr << endl << "**ERROR** in averaged: vector sizes do not match!" << endl << endl; + return 0.0; + } + + double sum_n = 0.0; + double result = 0.0; + for (unsigned int i=0; i + * + * \param fln file name + * \param param photo Meissner parameter set used. + * \param zz position vector (in nm) + * \param bb photo Meissner screening profile + * \param bbL London Meissner screening profile + */ +int dumpFile(const char *fln, const PhotoParam param, vector &zz, vector &bb, vector &bbL) +{ + ofstream fout(fln); + if (!fout.is_open()) { + cerr << endl << "**ERROR** couldn't open '" << fln << "' for writting." << endl << endl; + return 1; + } + + // write header + if (param.filmThickness == -1.0) + fout << "% semi-infinite" << endl; + else + fout << "% film with thickness: " << param.filmThickness << " (nm)" << endl; + fout << "% Parameters:" << endl; + fout << "% lamL = " << param.lamL << " (nm)" << endl; + fout << "% lamP = " << param.lamP << " (nm)" << endl; + fout << "% z0 = " << param.z0 << " (nm)" << endl; + fout << "% deadLayer = " << param.deadLayer << " (nm)" << endl; + if (zz.size() == 1) { + fout << "% --at =" << zz[0] << " (nm) used." << endl; + } else { + fout << "% --range =" << zz[0] << ", =" << zz[zz.size()-1] << " in (nm) used." << endl; + } + if (param.filmThickness == -1.0) { + fout << "% z B/Bext exp(-(z-deadLayer)/lamL) B/Bext-exp(-(z-deadLayer)/lamL)" << endl; + } else { + fout << "% z B/Bext cosh[((t-2*deadLayer)-z)/lamL]/cosh[(t-2*deadLayer)/lamL]=:f(z) B/Bext-f(z)" << endl; + } + for (int i=0; i + * + * \param fln file name + * \param param photo Meissner parameter set used. + * \param rgeFln rge-path-file-name vector + * \param zz averaged position vector + * \param bb averaged photo Meissner screening vector + * \param bbL London Meissner screening vector + */ +int dumpAvgFile(const char *fln, const PhotoParam param, const vector rgeFln, vector &zz, vector &bb, vector &bbL) +{ + ofstream fout(fln); + if (!fout.is_open()) { + cerr << endl << "**ERROR** couldn't open '" << fln << "' for writting." << endl << endl; + return 1; + } + + // write header + if (param.filmThickness == -1.0) + fout << "% semi-infinite" << endl; + else + fout << "% film with thickness: " << param.filmThickness << " (nm)" << endl; + fout << "% Parameters:" << endl; + fout << "% lamL = " << param.lamL << " (nm)" << endl; + fout << "% lamP = " << param.lamP << " (nm)" << endl; + fout << "% z0 = " << param.z0 << " (nm)" << endl; + fout << "% deadLayer = " << param.deadLayer << " (nm)" << endl; + fout << "% rge-file list:" << endl; + for (unsigned int i=0; i /Bext /Bext --- B: photo Meissner, BL: London Meissner" << endl; + for (unsigned int i=0; i irregular negative bessel function + * + * \param nu + * \param x + */ double InuMinus(double nu, double x) { return gsl_sf_bessel_Inu(nu,x) + M_2_PI * sin(nu*M_PI) * gsl_sf_bessel_Knu(nu,x); @@ -72,17 +569,28 @@ double InuMinus(double nu, double x) */ int main(int argc, char *argv[]) { - double lamL = -1.0; - double lamP = -1.0; - double z0 = -1.0; - double deadLayer = -1.0; - double filmThickness = -1.0; + PhotoParam param; + // init photo param + param.lamL = -1.0; + param.lamP = -1.0; + param.z0 = -1.0; + param.deadLayer = -1.0; + param.filmThickness = -1.0; + double zStart = -1.0; double zEnd = -1.0; char fln[1024]; + char rgeFln[1024]; + string prefix(""); + vector rgeFlnList; + vector< vector > rgeZ; + vector< vector > rgeN; + int jobTag=0; // 1=calc field without muon stopping distribution; 2=, for a given implantation energy (single rge-file); 3=, for list of rge-files strcpy(fln, ""); + strcpy(rgeFln, ""); + // handle input parameters if (argc <= 1) { photo_meissner_syntax(); return 0; @@ -101,50 +609,51 @@ int main(int argc, char *argv[]) int status=0; for (int i=1; i i+1) { status = sscanf(argv[i+1], "%lf", &zStart); if (status != 1) { - cout << endl << "**ERROR** --at =" << argv[i+1] << " is not a double." << endl; + cerr << endl << "**ERROR** --at =" << argv[i+1] << " is not a double." << endl; photo_meissner_syntax(); return 1; } } else { - cout << endl << "**ERROR** --at without found." << endl; + cerr << endl << "**ERROR** --at without found." << endl; photo_meissner_syntax(); return 1; } @@ -152,21 +661,22 @@ int main(int argc, char *argv[]) } if (!strcmp(argv[i], "--range")) { + jobTag = 2; if (argc > i+2) { status = sscanf(argv[i+1], "%lf", &zStart); if (status != 1) { - cout << endl << "**ERROR** --range : =" << argv[i+1] << " is not a double." << endl; + cerr << endl << "**ERROR** --range : =" << argv[i+1] << " is not a double." << endl; photo_meissner_syntax(); return 1; } status = sscanf(argv[i+2], "%lf", &zEnd); if (status != 1) { - cout << endl << "**ERROR** --range : =" << argv[i+2] << " is not a double." << endl; + cerr << endl << "**ERROR** --range : =" << argv[i+2] << " is not a double." << endl; photo_meissner_syntax(); return 1; } } else { - cout << endl << "**ERROR** --range without and/or found." << endl; + cerr << endl << "**ERROR** --range without and/or found." << endl; photo_meissner_syntax(); return 1; } @@ -177,165 +687,205 @@ int main(int argc, char *argv[]) if (argc > i+1) { strcpy(fln, argv[i+1]); } else { - cout << endl << "**ERROR** --dumpe without found." << endl; + cerr << endl << "**ERROR** --dump without found." << endl; photo_meissner_syntax(); return 1; } i += 1; } + + if (!strcmp(argv[i], "--file")) { + jobTag = 3; + if (argc > i+1) { + strcpy(rgeFln, argv[i+1]); + } else { + cerr << endl << "**ERROR** --file without found." << endl; + photo_meissner_syntax(); + return 1; + } + i += 1; + } + + if (!strcmp(argv[i], "--prefix")) { + jobTag = 4; + if (argc > i+1) { + prefix = argv[i+1]; + } else { + cerr << endl << "**ERROR** --prefix without found." << endl; + photo_meissner_syntax(); + return 1; + } + i += 1; + } + + if (!strcmp(argv[i], "--energies")) { + if (argc > i+1) { + // check if format is either a comma separated list, or of the form start:end:step + string str(argv[i+1]); + if (!generateRgeFln(prefix, str, rgeFlnList)) { + photo_meissner_syntax(); + return 1; + } + if (!readAllRgeFiles(rgeFlnList, rgeZ, rgeN)) + return 1; + } else { + cerr << endl << "**ERROR** --energies without found." << endl; + photo_meissner_syntax(); + return 1; + } + } } + // read single rge-file + if (strlen(rgeFln) > 0) { + vector zz, nn; + if (!readRgeFile(rgeFln, zz, nn)) { + return 1; + } + rgeZ.clear(); + rgeN.clear(); + rgeZ.resize(1); + rgeN.resize(1); + rgeZ[0] = zz; + rgeN[0] = nn; + } - if ((lamL == -1.0) || (lamP == -1.0) || (z0 == -1.0) || (deadLayer == -1.0) || (zStart == -1.0)) { - cout << endl << "**ERROR** not all needed parameters found. All the parameters need to be positive." << endl; + // check if all necessary parameters are present + if ((param.lamL == -1.0) || (param.lamP == -1.0) || (param.z0 == -1.0) || (param.deadLayer == -1.0) || ((zStart == -1.0) && (rgeZ.size() == 0))) { + cerr << endl << "**ERROR** not all needed parameters found. All the parameters need to be positive." << endl; photo_meissner_syntax(); return 1; } - if (lamL == 0.0) { - cout << endl << "**ERROR** lamL > 0.0 needed!" << endl << endl; + if (param.lamL == 0.0) { + cerr << endl << "**ERROR** lamL > 0.0 needed!" << endl << endl; return 1; } - if (lamP == 0.0) { - cout << endl << "**ERROR** lamP > 0.0 needed!" << endl << endl; + if (param.lamP == 0.0) { + cerr << endl << "**ERROR** lamP > 0.0 needed!" << endl << endl; return 1; } - if (z0 == 0.0) { - cout << endl << "**ERROR** z0 > 0.0 needed!" << endl << endl; + if (param.z0 == 0.0) { + cerr << endl << "**ERROR** z0 > 0.0 needed!" << endl << endl; return 1; } - if (filmThickness == 0.0) { - cout << endl << "**ERROR** filmThickness > 0.0 needed!" << endl << endl; + if (param.filmThickness == 0.0) { + cerr << endl << "**ERROR** filmThickness > 0.0 needed!" << endl << endl; return 1; } - // estimate dz - double dz = 0.0; - if (lamL < lamP) - dz = lamL / 500.0; - else - dz = lamP / 500.0; - - vector zz; - vector BB; - vector BBL; - double b, b1, bL; // b = B/Bext, bL = exp(-z/lamL) - if (filmThickness == -1.0) { // semi-infinite - double nuPlus = 2.0*z0/lamL; - double nuPhoto = 2.0*z0/lamP; - if (zEnd == -1.0) { // only a single point - if (zStart < deadLayer) { - b = 1.0; - bL = 1.0; - } else { - b1 = gsl_sf_bessel_Inu(nuPlus, nuPhoto); - assert(b1>0.0); - b = gsl_sf_bessel_Inu(nuPlus, nuPhoto * sqrt(exp(-(zStart-deadLayer)/z0)))/b1; - bL = exp(-(zStart-deadLayer)/lamL); - } - - if (!strcmp(fln, "")) // only to stdout if no dump present - cout << endl << "result: zPos=" << zStart << ", B/Bext=" << b << ", exp(-[zPos-deadLayer]/lamL)=" << bL << endl << endl; - // keep z, b - zz.push_back(zStart); - BB.push_back(b); - BBL.push_back(bL); - } else { // a range - double z = zStart; - b1 = gsl_sf_bessel_Inu(nuPlus, nuPhoto); - assert(b1>0.0); - do { - if (z < deadLayer) { - b = 1.0; + // do the calculation depending on the job request + double b=-1.0, bL = -1.0; + vector zz, bb, bbL, nn; + vector zMean, BMean, BLMean; + switch (jobTag) { + case 1: // single field value + b = calcSingleField(zStart, param); + if (param.filmThickness == -1.0) { // semi-infinite + if (zStart < param.deadLayer) bL = 1.0; - } else { - b = gsl_sf_bessel_Inu(nuPlus, nuPhoto * sqrt(exp(-(z-deadLayer)/z0)))/b1; - bL = exp(-(z-deadLayer)/lamL); - } - if (!strcmp(fln, "")) // only to stdout if no dump present - cout << z << ", " << b << ", " << bL << ", " << b-bL << endl; - zz.push_back(z); - BB.push_back(b); - BBL.push_back(bL); - z += dz; - } while (z <= zEnd); - } - } else { // film - double nuPlus = 2.0*z0/lamL; - double nuPhoto = 2.0*z0/lamP; - double beta = exp(-(filmThickness/2.0-deadLayer)/z0); - double NN = InuMinus(nuPlus, nuPhoto*beta)*gsl_sf_bessel_Inu(nuPlus, nuPhoto) - - InuMinus(nuPlus, nuPhoto)*gsl_sf_bessel_Inu(nuPlus, nuPhoto*beta); - assert(NN>0); - if (zEnd == -1.0) { // only a single point - if ((zStart < deadLayer) || (zStart > filmThickness-deadLayer)) { - b = 1.0; - bL = 1.0; + else + bL = exp(-(zStart-param.deadLayer)/param.lamL); } else { - b = (InuMinus(nuPlus, nuPhoto*exp(-(zStart-deadLayer)/(2.0*z0)))*(gsl_sf_bessel_Inu(nuPlus, nuPhoto)-gsl_sf_bessel_Inu(nuPlus, nuPhoto*beta))- - gsl_sf_bessel_Inu(nuPlus, nuPhoto*exp(-(zStart-deadLayer)/(2.0*z0)))*(InuMinus(nuPlus, nuPhoto)-InuMinus(nuPlus, nuPhoto*beta)))/NN; - bL = cosh((filmThickness/2.0-zStart)/lamL)/cosh((filmThickness/2.0-deadLayer)/lamL); - } - - if (!strcmp(fln, "")) // only to stdout if no dump present - cout << endl << "result: zPos=" << zStart << ", B/Bext=" << b << ", cosh[((t-deadLayer)-z)/lamL]/cosh[(t-deadLayer)/lamL]=:f(z)=" << bL << ", B/Bext-f(z)=" << b-bL << endl << endl; - // keep z, b - zz.push_back(zStart); - BB.push_back(b); - BBL.push_back(bL); - } else { // a range - double z = zStart; - do { - if ((z < deadLayer) || (z > filmThickness-deadLayer)) { - b = 1.0; + if ((zStart < param.deadLayer) || (zStart > param.filmThickness - param.deadLayer)) bL = 1.0; - } else { - b = (InuMinus(nuPlus, nuPhoto*exp(-(z-deadLayer)/(2.0*z0)))*(gsl_sf_bessel_Inu(nuPlus, nuPhoto)-gsl_sf_bessel_Inu(nuPlus, nuPhoto*beta))- - gsl_sf_bessel_Inu(nuPlus, nuPhoto*exp(-(z-deadLayer)/(2.0*z0)))*(InuMinus(nuPlus, nuPhoto)-InuMinus(nuPlus, nuPhoto*beta)))/NN; - bL = cosh((filmThickness/2.0-z)/lamL)/cosh((filmThickness/2.0-deadLayer)/lamL); + else + bL = cosh((param.filmThickness/2.0-zStart)/param.lamL)/cosh((param.filmThickness/2.0-param.deadLayer)/param.lamL); + } + if (!strcmp(fln, "")) { // only to stdout if no dump present + if (param.filmThickness == -1.0) // semi-infinite + cout << endl << "result: zPos=" << zStart << ", B/Bext=" << b << ", exp(-[zPos-deadLayer]/lamL)=:f(z)=" << bL << ", B/Bext-f(z)=" << b-bL << endl << endl; + else // film + cout << endl << "result: zPos=" << zStart << ", B/Bext=" << b << ", cosh[((t-deadLayer)-z)/lamL]/cosh[(t-deadLayer)/lamL]=:f(z)=" << bL << ", B/Bext-f(z)=" << b-bL << endl << endl; + } else { + zz.push_back(zStart); + bb.push_back(b); + bbL.push_back(bL); + dumpFile(fln, param, zz, bb, bbL); + } + break; + case 2: // range of field values + calcFields(zStart, zEnd, param, zz, bb, bbL); + if (!strcmp(fln, "")) { // only to stdout if no dump present + // write header + if (param.filmThickness == -1.0) + cout << "% semi-infinite" << endl; + else + cout << "% film with thickness: " << param.filmThickness << " (nm)" << endl; + cout << "% Parameters:" << endl; + cout << "% lamL = " << param.lamL << " (nm)" << endl; + cout << "% lamP = " << param.lamP << " (nm)" << endl; + cout << "% z0 = " << param.z0 << " (nm)" << endl; + cout << "% deadLayer = " << param.deadLayer << " (nm)" << endl; + cout << "% --range =" << zz[0] << ", =" << zz[zz.size()-1] << " in (nm) used." << endl; + if (param.filmThickness == -1.0) + cout << "% z B/Bext exp(-(z-deadLayer)/lamL) B/Bext-exp(-(z-deadLayer)/lamL)" << endl; + else + cout << "% z B/Bext cosh[((t-2*deadLayer)-z)/lamL]/cosh[(t-2*deadLayer)/lamL]=:f(z) B/Bext-f(z)" << endl; + for (int i=0; i=" << zStart << " (nm) used." << endl; - } else { - fout << "% --range =" << zStart << ", =" << zEnd << " in (nm) used." << endl; - } - if (filmThickness == -1.0) - fout << "% z B/Bext exp(-(z-deadLayer)/lamL) B/Bext-exp(-(z-deadLayer)/lamL)" << endl; - else - fout << "% z B/Bext cosh[((t-2*deadLayer)-z)/lamL]/cosh[(t-2*deadLayer)/lamL]=:f(z) B/Bext-f(z)" << endl; - for (int i=0; i, for a single rge-file + // define a proper z-range + zStart = 0.0; + if (param.filmThickness == -1.0) + zEnd = 10.0*param.lamL; + else + zEnd = param.filmThickness; + // calculate the field for the given parameters + calcFields(zStart, zEnd, param, zz, bb, bbL); + // generate a matched n(z) vector + nn.clear(); + matched_nz(rgeZ[0], rgeN[0], zz, nn); + zMean.push_back(averaged(zz, nn)); + BMean.push_back(averaged(zz, bb)); + BLMean.push_back(averaged(zz, bbL)); + rgeFlnList.push_back(rgeFln); + if (!strcmp(fln, "")) + cout << "=" << zMean[0] << " (nm), /Bext=" << BMean[0] << ", /Bext=" << BLMean[0] << " for rge-file: '" << rgeFln << "'." << endl; + else + dumpAvgFile(fln, param, rgeFlnList, zMean, BMean, BLMean); + break; + case 4: // , for a list of rge-files + // define a proper z-range + zStart = 0.0; + if (param.filmThickness == -1.0) + zEnd = 10.0*param.lamL; + else + zEnd = param.filmThickness; + // calculate the field for the given parameters + calcFields(zStart, zEnd, param, zz, bb, bbL); + for (unsigned int i=0; i (nm), /Bext, /Bext" << endl; + for (unsigned int i=0; i