Start writing the general inf. lat. utility

Bulk copied from 2dlattice.c + use gengetopts


Former-commit-id: 27a21038f292dc10bfbb36f3dfd509f2d0871240
This commit is contained in:
Marek Nečada 2019-06-18 13:36:57 +03:00
parent c4e617a5b1
commit 96832355c7
5 changed files with 1518 additions and 2 deletions

View File

@ -25,4 +25,4 @@ use_c99()
# ) # )
#target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) #target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
add_executable(EwaldTransOps transop-ewald.c transop_ewald_cmdline.c) add_executable(EwaldTransOps transop_ewald.c transop_ewald_cmdline.c)

View File

@ -14,7 +14,11 @@ option "output" o "Output file"
default="" default=""
optional optional
option "normalisation" n "VSWF normalisation convention" option "error-estimate-output" E "Path to the output with error estimates"
string
optional
option "normalisation" N "VSWF normalisation convention"
values="Power","None","SH" enum values="Power","None","SH" enum
default="Power" default="Power"
@ -27,12 +31,47 @@ option "Ewald-parameter" e "The value of Ewald parameter η"
double double
optional optional
option "frequency-unit" u "Specifies the frequency unit is used for inputs."
values="eV","scuff"
default="scuff"
option "lMax" L "Maximum spherical multipole order to which the translation operator elements are calculated"
int
required
option "eta" n "Medium refractive index"
double
required
option "particle" p "Specify the x and y coordinates of a single particle; If not specified, one particle per unit cell is assumed."
optional
multiple
defmode "k_omega_points" modedesc="Specifying each (ω, k) pair separately." defmode "k_omega_points" modedesc="Specifying each (ω, k) pair separately."
defmode "k_omega_meshgrid" modedesc="Specifying lists of ω and k, from which all possible pairs are generated." defmode "k_omega_meshgrid" modedesc="Specifying lists of ω and k, from which all possible pairs are generated."
modeoption "pointfile" T "Path to a file containing frequency, k_x, k_y triples\
(separated by white spaces). If not specified, read them from stdin."
mode="k_omega_points" multiple default="-"
modeoption "point" t "Specifies a frequency, k_x, k_y triple, separated by commas."
mode="k_omega_points" multiple optional
modeoption "omegafile" F "Path to a file containing a list of frequencies\
separated by whitespaces."
mode="k_omega_meshgrid" multiple optional
modeoption "omega" f "Specifies frequency (or multiple frequencies separated by semicolons) on the command line."
mode="k_omega_meshgrid" multiple optional
modeoption "kfile" K "Path to a file containing a list of k_x, k_y pairs."
mode="k_omega_meshgrid" multiple optional
default="-"
modeoption "k" k "Specifies pair(s) of k_x, k_y values"
mode="k_omega_meshgrid" multiple optional
#option <long> <short> "<desc>" #option <long> <short> "<desc>"
# {details="<detailed description>"} # {details="<detailed description>"}
# {argtype} {typestr="<type descr>"} # {argtype} {typestr="<type descr>"}

View File

@ -0,0 +1,160 @@
// c99 -o ew_gen_kin -Wall -I ../.. -I ../../amos/ -O2 -ggdb -DQPMS_VECTORS_NICE_TRANSFORMATIONS -DLATTICESUMS32 2dlattice_ewald.c ../translations.c ../ewald.c ../ewaldsf.c ../gaunt.c ../lattices2d.c ../latticegens.c ../bessel.c -lgsl -lm -lblas ../../amos/libamos.a -lgfortran ../error.c
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "transop_ewald_cmdline.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#define LATTICESUMS32
#include <qpms/translations.h>
#include <qpms/lattices.h>
#include <gsl/gsl_const_mksa.h>
// Command line order:
// outfile b1.x b1.y b2.x b2.y lMax scuffomega refindex npart part0.x part0.y [part1.x part1.y [...]]
//
// Standard input (per line):
// k.x k.y
//
// Output data format (line):
//
#define MAXKCOUNT 200 // 200 // serves as klist default buffer size
//#define KMINCOEFF 0.783 //0.9783 // 0.783 // not used if KSTDIN defined
//#define KMAXCOEFF 1.217 //1.0217 // 1.217 // not used if KSTDIN defined
#define KLAYERS 20
#define RLAYERS 20
const double s3 = 1.732050807568877293527446341505872366942805253810380628055;
//const qpms_y_t lMax = 3;
//const double REFINDEX = 1.52;
static const double SCUFF_OMEGAUNIT = 3e14;
static const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
static const double eV = GSL_CONST_MKSA_ELECTRON_CHARGE;
static const double c0 = GSL_CONST_MKSA_SPEED_OF_LIGHT;
int main (int argc, char **argv) {
struct gengetopt_args_info args_info;
int retval = cmdline_parser(argc, argv, *argc_info);
if (retval) return retval;
char *outfile = argv[1];
char *errfile = NULL; // Filename for the error estimate output; NOT USED
cart2_t b1 = {strtod(argv[2], NULL), strtod(argv[3], NULL)},
b2 = {strtod(argv[4], NULL), strtod(argv[5], NULL)};
const qpms_l_t lMax = strtol(argv[6], NULL, 10); assert(lMax>0);
const double scuffomega = strtod(argv[7], NULL);
const double refindex = strtod(argv[8], NULL);
const int npart = strtol(argv[9], NULL, 10); assert(npart>0);
assert(argc >= 2*npart + 10);
assert(npart > 0);
cart2_t part_positions[npart];
for(int p = 0; p < npart; ++p) {
part_positions[p].x = strtod(argv[10+2*p], NULL);
part_positions[p].y = strtod(argv[10+2*p+1], NULL);
}
//#ifdef KSTDIN
size_t kcount = 0;
size_t klist_capacity = MAXKCOUNT;
cart2_t *klist = malloc(sizeof(cart2_t) * klist_capacity);
while (scanf("%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) {
++kcount;
if(kcount >= klist_capacity) {
klist_capacity *= 2;
klist = realloc(klist, sizeof(cart2_t) * klist_capacity);
if (klist == NULL) abort();
}
}
//#else
#if 0
cart2_t klist[MAXKCOUNT];
int kcount = MAXKCOUNT;
for (int i = 0; i < kcount; ++i) { // TODO this should depend on orientation...
klist[i].x = 0;
klist[i].y = (4.* M_PI / 3. / LATTICE_A) * (KMINCOEFF + (KMAXCOEFF-KMINCOEFF)/kcount*i);
}
#endif
const double unitcell_area = l2d_unitcell_area(b1, b2);
l2d_reduceBasis(b1, b2, &b1, &b2);
// TODO more clever way of determining the cutoff
const double a = sqrt(unitcell_area); // N.B. different meaning than before
const double maxR = 25 * a;
const double maxK = 25 * 2*M_PI/a;
qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_POWER_CS); // vai POWER_CS?
FILE *out = fopen(outfile, "w");
FILE *err = NULL;
if (errfile)
err = fopen(errfile, "w");
{
const double omega = scuffomega * SCUFF_OMEGAUNIT;
const double EeV = omega * hbar / eV;
const double k0_vac = omega / c0;
const double k0_eff = k0_vac * refindex;
const double eta = 5.224/a; // FIXME quite arbitrary, but this one should work
// indices : destpart (A/B-particle), srcpart (A/B-particle), coeff type (A/B- type), desty, srcy
complex double W[npart][npart][2][c->nelem][c->nelem];
double Werr[npart][npart][npart][c->nelem][c->nelem];
for (size_t ki = 0; ki < kcount; ++ki) {
cart2_t beta = klist[ki];
memset(W, 0, sizeof(W));
if(err)
memset(Werr, 0, sizeof(Werr));
const ptrdiff_t deststride = &(W[0][0][0][1][0]) - &(W[0][0][0][0][0]);
const ptrdiff_t srcstride = &(W[0][0][0][0][1]) - &(W[0][0][0][0][0]);
assert (srcstride == 1 && deststride == c->nelem);
for (size_t ps = 0; ps < npart; ++ps) for (size_t pd = 0; pd < npart; ++pd)
// TODO optimize (calculate only once for each particle shift; especially if pd == ps)
qpms_trans_calculator_get_AB_arrays_e32(c,
&(W[pd][ps][0][0][0]), err ? &(Werr[pd][ps][0][0][0]) : NULL, // Adest, Aerr,
&(W[pd][ps][1][0][0]), err ? &(Werr[pd][ps][1][0][0]) : NULL, // Bdest, Berr,
deststride, srcstride,
eta, k0_eff, b1, b2,
beta,
cart2_substract(part_positions[pd], part_positions[ps]), // CHECKSIGN
maxR, maxK
);
// TODO CHECK B<-A vs. A<-B relation
fprintf(out, "%.16g\t%.16g\t%.16g\t%.16g\t%.16g\t",
scuffomega, EeV, k0_eff, beta.x, beta.y);
if(err) fprintf(err, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
scuffomega, EeV, k0_eff, beta.x, beta.y);
size_t totalelems = sizeof(W) / sizeof(complex double);
for (size_t i = 0; i < totalelems; ++i) {
complex double w = ((complex double *)W)[i];
fprintf(out, "%.16g\t%.16g\t", creal(w), cimag(w));
if (err)
fprintf(err, "%.3g\t", ((double *)Werr)[i]);
}
fputc('\n', out);
if(err) fputc('\n', err);
}
}
fclose(out);
if(err) fclose(err);
//#ifdef KSTDIN
free(klist);
//#endif
qpms_trans_calculator_free(c);
}

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,249 @@
/** @file transop_ewald_cmdline.h
* @brief The header file for the command line option parser
* generated by GNU Gengetopt version 2.22.6
* http://www.gnu.org/software/gengetopt.
* DO NOT modify this file, since it can be overwritten
* @author GNU Gengetopt by Lorenzo Bettini */
#ifndef TRANSOP_EWALD_CMDLINE_H
#define TRANSOP_EWALD_CMDLINE_H
/* If we use autoconf. */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h> /* for FILE */
#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */
#ifndef CMDLINE_PARSER_PACKAGE
/** @brief the program name (used for printing errors) */
#define CMDLINE_PARSER_PACKAGE "qpms-translations-ewald32"
#endif
#ifndef CMDLINE_PARSER_PACKAGE_NAME
/** @brief the complete program name (used for help and version) */
#define CMDLINE_PARSER_PACKAGE_NAME "qpms-translations-ewald32"
#endif
#ifndef CMDLINE_PARSER_VERSION
/** @brief the program version */
#define CMDLINE_PARSER_VERSION "dev"
#endif
enum enum_normalisation { normalisation__NULL = -1, normalisation_arg_Power = 0, normalisation_arg_None, normalisation_arg_SH };
/** @brief Where the command line options are stored */
struct gengetopt_args_info
{
const char *help_help; /**< @brief Print help and exit help description. */
const char *detailed_help_help; /**< @brief Print help, including all details and hidden options, and exit help description. */
const char *version_help; /**< @brief Print version and exit help description. */
char * output_arg; /**< @brief Output file (default=''). */
char * output_orig; /**< @brief Output file original value given at command line. */
const char *output_help; /**< @brief Output file help description. */
char * error_estimate_output_arg; /**< @brief Path to the output with error estimates. */
char * error_estimate_output_orig; /**< @brief Path to the output with error estimates original value given at command line. */
const char *error_estimate_output_help; /**< @brief Path to the output with error estimates help description. */
enum enum_normalisation normalisation_arg; /**< @brief VSWF normalisation convention (default='Power'). */
char * normalisation_orig; /**< @brief VSWF normalisation convention original value given at command line. */
const char *normalisation_help; /**< @brief VSWF normalisation convention help description. */
int csphase_arg; /**< @brief Whether the Condon-Shortley phase is included in VSWF definition (-1) or not (+1) (default='-1'). */
char * csphase_orig; /**< @brief Whether the Condon-Shortley phase is included in VSWF definition (-1) or not (+1) original value given at command line. */
const char *csphase_help; /**< @brief Whether the Condon-Shortley phase is included in VSWF definition (-1) or not (+1) help description. */
double Ewald_parameter_arg; /**< @brief The value of Ewald parameter η. */
char * Ewald_parameter_orig; /**< @brief The value of Ewald parameter η original value given at command line. */
const char *Ewald_parameter_help; /**< @brief The value of Ewald parameter η help description. */
char * frequency_unit_arg; /**< @brief Specifies the frequency unit is used for inputs. (default='scuff'). */
char * frequency_unit_orig; /**< @brief Specifies the frequency unit is used for inputs. original value given at command line. */
const char *frequency_unit_help; /**< @brief Specifies the frequency unit is used for inputs. help description. */
int lMax_arg; /**< @brief Maximum spherical multipole order to which the translation operator elements are calculated. */
char * lMax_orig; /**< @brief Maximum spherical multipole order to which the translation operator elements are calculated original value given at command line. */
const char *lMax_help; /**< @brief Maximum spherical multipole order to which the translation operator elements are calculated help description. */
double eta_arg; /**< @brief Medium refractive index. */
char * eta_orig; /**< @brief Medium refractive index original value given at command line. */
const char *eta_help; /**< @brief Medium refractive index help description. */
unsigned int particle_min; /**< @brief Specify the x and y coordinates of a single particle; If not specified, one particle per unit cell is assumed.'s minimum occurreces */
unsigned int particle_max; /**< @brief Specify the x and y coordinates of a single particle; If not specified, one particle per unit cell is assumed.'s maximum occurreces */
const char *particle_help; /**< @brief Specify the x and y coordinates of a single particle; If not specified, one particle per unit cell is assumed. help description. */
unsigned int pointfile_min; /**< @brief Path to a file containing frequency, k_x, k_y triples(separated by white spaces). If not specified, read them from stdin.'s minimum occurreces */
unsigned int pointfile_max; /**< @brief Path to a file containing frequency, k_x, k_y triples(separated by white spaces). If not specified, read them from stdin.'s maximum occurreces */
const char *pointfile_help; /**< @brief Path to a file containing frequency, k_x, k_y triples(separated by white spaces). If not specified, read them from stdin. help description. */
unsigned int point_min; /**< @brief Specifies a frequency, k_x, k_y triple, separated by commas.'s minimum occurreces */
unsigned int point_max; /**< @brief Specifies a frequency, k_x, k_y triple, separated by commas.'s maximum occurreces */
const char *point_help; /**< @brief Specifies a frequency, k_x, k_y triple, separated by commas. help description. */
unsigned int omegafile_min; /**< @brief Path to a file containing a list of frequenciesseparated by whitespaces.'s minimum occurreces */
unsigned int omegafile_max; /**< @brief Path to a file containing a list of frequenciesseparated by whitespaces.'s maximum occurreces */
const char *omegafile_help; /**< @brief Path to a file containing a list of frequenciesseparated by whitespaces. help description. */
unsigned int omega_min; /**< @brief Specifies frequency (or multiple frequencies separated by semicolons) on the command line.'s minimum occurreces */
unsigned int omega_max; /**< @brief Specifies frequency (or multiple frequencies separated by semicolons) on the command line.'s maximum occurreces */
const char *omega_help; /**< @brief Specifies frequency (or multiple frequencies separated by semicolons) on the command line. help description. */
unsigned int kfile_min; /**< @brief Path to a file containing a list of k_x, k_y pairs.'s minimum occurreces */
unsigned int kfile_max; /**< @brief Path to a file containing a list of k_x, k_y pairs.'s maximum occurreces */
const char *kfile_help; /**< @brief Path to a file containing a list of k_x, k_y pairs. help description. */
unsigned int k_min; /**< @brief Specifies pair(s) of k_x, k_y values's minimum occurreces */
unsigned int k_max; /**< @brief Specifies pair(s) of k_x, k_y values's maximum occurreces */
const char *k_help; /**< @brief Specifies pair(s) of k_x, k_y values help description. */
unsigned int help_given ; /**< @brief Whether help was given. */
unsigned int detailed_help_given ; /**< @brief Whether detailed-help was given. */
unsigned int version_given ; /**< @brief Whether version was given. */
unsigned int output_given ; /**< @brief Whether output was given. */
unsigned int error_estimate_output_given ; /**< @brief Whether error-estimate-output was given. */
unsigned int normalisation_given ; /**< @brief Whether normalisation was given. */
unsigned int csphase_given ; /**< @brief Whether csphase was given. */
unsigned int Ewald_parameter_given ; /**< @brief Whether Ewald-parameter was given. */
unsigned int frequency_unit_given ; /**< @brief Whether frequency-unit was given. */
unsigned int lMax_given ; /**< @brief Whether lMax was given. */
unsigned int eta_given ; /**< @brief Whether eta was given. */
unsigned int particle_given ; /**< @brief Whether particle was given. */
unsigned int pointfile_given ; /**< @brief Whether pointfile was given. */
unsigned int point_given ; /**< @brief Whether point was given. */
unsigned int omegafile_given ; /**< @brief Whether omegafile was given. */
unsigned int omega_given ; /**< @brief Whether omega was given. */
unsigned int kfile_given ; /**< @brief Whether kfile was given. */
unsigned int k_given ; /**< @brief Whether k was given. */
int k_omega_meshgrid_mode_counter; /**< @brief Counter for mode k_omega_meshgrid */
int k_omega_points_mode_counter; /**< @brief Counter for mode k_omega_points */
} ;
/** @brief The additional parameters to pass to parser functions */
struct cmdline_parser_params
{
int override; /**< @brief whether to override possibly already present options (default 0) */
int initialize; /**< @brief whether to initialize the option structure gengetopt_args_info (default 1) */
int check_required; /**< @brief whether to check that all required options were provided (default 1) */
int check_ambiguity; /**< @brief whether to check for options already specified in the option structure gengetopt_args_info (default 0) */
int print_errors; /**< @brief whether getopt_long should print an error message for a bad option (default 1) */
} ;
/** @brief the purpose string of the program */
extern const char *gengetopt_args_info_purpose;
/** @brief the usage string of the program */
extern const char *gengetopt_args_info_usage;
/** @brief the description string of the program */
extern const char *gengetopt_args_info_description;
/** @brief all the lines making the help output */
extern const char *gengetopt_args_info_help[];
/** @brief all the lines making the detailed help output (including hidden options and details) */
extern const char *gengetopt_args_info_detailed_help[];
/**
* The command line parser
* @param argc the number of command line options
* @param argv the command line options
* @param args_info the structure where option information will be stored
* @return 0 if everything went fine, NON 0 if an error took place
*/
int cmdline_parser (int argc, char **argv,
struct gengetopt_args_info *args_info);
/**
* The command line parser (version with additional parameters - deprecated)
* @param argc the number of command line options
* @param argv the command line options
* @param args_info the structure where option information will be stored
* @param override whether to override possibly already present options
* @param initialize whether to initialize the option structure my_args_info
* @param check_required whether to check that all required options were provided
* @return 0 if everything went fine, NON 0 if an error took place
* @deprecated use cmdline_parser_ext() instead
*/
int cmdline_parser2 (int argc, char **argv,
struct gengetopt_args_info *args_info,
int override, int initialize, int check_required);
/**
* The command line parser (version with additional parameters)
* @param argc the number of command line options
* @param argv the command line options
* @param args_info the structure where option information will be stored
* @param params additional parameters for the parser
* @return 0 if everything went fine, NON 0 if an error took place
*/
int cmdline_parser_ext (int argc, char **argv,
struct gengetopt_args_info *args_info,
struct cmdline_parser_params *params);
/**
* Save the contents of the option struct into an already open FILE stream.
* @param outfile the stream where to dump options
* @param args_info the option struct to dump
* @return 0 if everything went fine, NON 0 if an error took place
*/
int cmdline_parser_dump(FILE *outfile,
struct gengetopt_args_info *args_info);
/**
* Save the contents of the option struct into a (text) file.
* This file can be read by the config file parser (if generated by gengetopt)
* @param filename the file where to save
* @param args_info the option struct to save
* @return 0 if everything went fine, NON 0 if an error took place
*/
int cmdline_parser_file_save(const char *filename,
struct gengetopt_args_info *args_info);
/**
* Print the help
*/
void cmdline_parser_print_help(void);
/**
* Print the detailed help (including hidden options and details)
*/
void cmdline_parser_print_detailed_help(void);
/**
* Print the version
*/
void cmdline_parser_print_version(void);
/**
* Initializes all the fields a cmdline_parser_params structure
* to their default values
* @param params the structure to initialize
*/
void cmdline_parser_params_init(struct cmdline_parser_params *params);
/**
* Allocates dynamically a cmdline_parser_params structure and initializes
* all its fields to their default values
* @return the created and initialized cmdline_parser_params structure
*/
struct cmdline_parser_params *cmdline_parser_params_create(void);
/**
* Initializes the passed gengetopt_args_info structure's fields
* (also set default values for options that have a default)
* @param args_info the structure to initialize
*/
void cmdline_parser_init (struct gengetopt_args_info *args_info);
/**
* Deallocates the string fields of the gengetopt_args_info structure
* (but does not deallocate the structure itself)
* @param args_info the structure to deallocate
*/
void cmdline_parser_free (struct gengetopt_args_info *args_info);
/**
* Checks that all the required options were specified
* @param args_info the structure to check
* @param prog_name the name of the program that will be used to print
* possible errors
* @return
*/
int cmdline_parser_required (struct gengetopt_args_info *args_info,
const char *prog_name);
extern const char *cmdline_parser_normalisation_values[]; /**< @brief Possible values for normalisation. */
extern const char *cmdline_parser_csphase_values[]; /**< @brief Possible values for csphase. */
extern const char *cmdline_parser_frequency_unit_values[]; /**< @brief Possible values for frequency-unit. */
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif /* TRANSOP_EWALD_CMDLINE_H */