apps/transop-ewald/transop_ewald.c

   1 // 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
   2
   3 #ifdef HAVE_CONFIG_H
   4 #include <config.h>
   5 #endif
   6
   7 #include "transop_ewald_cmdline.h"
   8
   9 #include <stdio.h>
  10 #include <math.h>
  11 #include <string.h>
  12 #include <errno.h>
  13 #define LATTICESUMS32
  14 #include <qpms/translations.h>
  15 #include <qpms/lattices.h>
  16 #include <qpms/qpms_error.h>
  17 #include <gsl/gsl_const_mksa.h>
  18 #include <qpms/parsing.h>
  19
  20
  21 // Command line args parsing progress:
  22 // output
  23 // base-vector DONE 2D
  24 // error-estimate-output
  25 // normalisation
  26 // csphase
  27 // Ewald-parameter
  28 // frequency-unit
  29 // lMax DONE
  30 // refractive-index DONE
  31 // particle DONE
  32 // pointfile
  33 // point
  34 // omegafile DONE, TODO unit conversion
  35 // omega DONE, TODO unit conversion
  36 // kfile DONE 2D
  37 // k DONE 2D
  38
  39 #define MAXKCOUNT 200 // 200 // serves as klist default buffer size
  40 //#define KMINCOEFF 0.783 //0.9783 // 0.783 // not used if KSTDIN defined
  41 //#define KMAXCOEFF 1.217 //1.0217 // 1.217 // not used if KSTDIN defined
  42 #define KLAYERS 20
  43 #define RLAYERS 20
  44
  45 const double s3 = 1.732050807568877293527446341505872366942805253810380628055;
  46
  47 //const qpms_y_t lMax = 3;
  48 //const double REFINDEX = 1.52;
  49 static const double SCUFF_OMEGAUNIT = 3e14;
  50 static const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
  51 static const double eV = GSL_CONST_MKSA_ELECTRON_CHARGE;
  52 static const double c0 = GSL_CONST_MKSA_SPEED_OF_LIGHT;
  53
  54 int main (int argc, char **argv) {
  55   struct gengetopt_args_info args_info;
  56
  57   int retval = cmdline_parser(argc, argv, *args_info);
  58   if (retval) return retval;
  59
  60   // Parse lattice vectors
  61   const int latdim = args_info.base_vector_given;
  62   QPMS_ENSURE(latdim == 2,
  63       "Sorry, only 2d lattices supported, but %d basis vectors were given\n",
  64       latdim);
  65   cart2_t b[latdim];
  66   for (int i = 0; i < latdim; ++i) {
  67     const int gotnumbers = qpms_parse_ndoubles(
  68           (*double) &(b[i].x), latdim,
  69           args_info.base_vector_arg[i]);
  70     QPMS_ENSURE(latdim == gotnumbers,
  71         "%d. base vector contained %d numbers, expected %d\n",
  72         i, gotnumbers, latdim);
  73   }
  74
  75   // N.B. this is 2D specific, TODO generalize when Nd sum supported
  76   const double unitcell_area = l2d_unitcell_area(b[0], b[1]);
  77   l2d_reduceBasis(b[0], b[1], b, b+1);
  78
  79   const qpms_l_t lMax = args_info.lMax_arg;
  80   QPMS_ENSURE(lMax > 0, "invalid value of lMax: %d", (int)lMax);
  81
  82   const double refindex = args_info.refractive_index_arg;
  83
  84   // Parse all particle positions
  85   const int npart = args_info.particle_given;
  86   if(!npart) ++npart;
  87   cart2_t part_positions[npart];
  88   if(!args_info.particle_given)
  89     part_positions[0].x = part_positions[0].y = 0;
  90   else for (int i = 0; i < npart; ++i) {
  91     const int gotnumbers = qpms_parse_ndoubles(
  92         (*double) &(part_positions[i].x), latdim,
  93         args_info.particle_arg[i]);
  94     QPMS_ENSURE(latdim == gotnumbers,
  95         "%d. particle position contained %d coordinates, expected %d\n",
  96         i, gotnumbers, latdim);
  97   }
  98
  99   QPMS_ENSURE(!args_info.k_omega_meshgrid_mode_counter !=
 100       !args_info.k_omega_points_mode_counter,
 101       "THIS IS A BUG. Only one mode ((k, ω) tuples, or k, ω lists) allowed.");
 102   // ===================== k, ω grid mode =====================
 103   if (args_info.k_omega_meshgrid_mode_counter) {
 104     size_t omegacount = 0;
 105     double *omegalist = NULL;
 106     for (int i = 0; i < args_info.omega_given; ++i) // freqs from command line
 107       omegacount = qpms_parse_doubles(&omegalist, omegacount,
 108           args_info.omega_arg[i]);
 109     for (int i = 0; i < args_info.omegafile_given; ++i) // freqs from file
 110       omegacount = qpms_parse_doubles_fromfile(&omegalist, omegacount,
 111           args_info.omegafile_arg[i]);
 112
 113     size_t kc_count = 0;
 114     double *kclist = NULL;
 115     for (int i = 0; i < args_info.k_given; ++i) {// ks from command line
 116       kc_count = qpms_parse_doubles(&kclist, kc_count, args_info.k_arg[i]);
 117       QPMS_ENSURE(0 == kc_count % latdim,
 118           "Provided number of k components (cum. %zd) not compatible with the "
 119           "lattice dimension (%d): %s", kc_count, latdim, args_info.k_arg[i]);
 120     }
 121     for (int i = 0; i < args_info.kfile_given; ++i) {//ks from file
 122       kc_count = qpms_parse_doubles_fromfile(&kclist, kc_count,
 123           args_info.kfile_arg[i]);
 124       QPMS_ENSURE(0 == kc_count % latdim,
 125           "Provided number of k components (cum. %zd) not compatible with the "
 126           "lattice dimension (%d) in file %s", kc_count, latdim,
 127           args_info.kfile_arg[i]);
 128     }
 129     // 2D specific, TODO generalize when Nd supported
 130     cart2_t klist[kc_count/2];
 131     for (size_t i = 0; i < kc_count/2; ++i)
 132       klist[i] = {kclist[2*i], kclist[2*i+1]};
 133     free(kclist);
 134
 135     TODO;
 136
 137   } else if (args_info.k_omega_points_mode_counter) { // explic. point mode
 138     TODO;
 139   }
 140
 141
 142   const double scuffomega = strtod(argv[7], NULL);
 143 //#ifdef KSTDIN
 144   size_t kcount = 0;
 145   size_t klist_capacity = MAXKCOUNT;
 146   cart2_t *klist = malloc(sizeof(cart2_t) * klist_capacity);
 147   while (scanf("%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) {
 148     ++kcount;
 149     if(kcount >= klist_capacity) {
 150       klist_capacity *= 2;
 151       klist = realloc(klist, sizeof(cart2_t) * klist_capacity);
 152       if (klist == NULL) abort();
 153     }
 154   }
 155 //#else
 156 #if 0
 157   cart2_t klist[MAXKCOUNT];
 158   int kcount = MAXKCOUNT;
 159   for (int i = 0; i <  kcount; ++i) { // TODO this should depend on orientation...
 160     klist[i].x = 0;
 161     klist[i].y = (4.* M_PI / 3. / LATTICE_A) * (KMINCOEFF + (KMAXCOEFF-KMINCOEFF)/kcount*i);
 162   }
 163 #endif
 164
 165
 166   // TODO more clever way of determining the cutoff
 167   const double a = sqrt(unitcell_area); // N.B. different meaning than before
 168   const double maxR = 25 * a;
 169   const double maxK = 25 * 2*M_PI/a;
 170
 171   qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_POWER_CS); // vai POWER_CS?
 172
 173   FILE *out, *ferr = NULL;
 174   if (args_info.error_estimate_output_given) {
 175     if (!strcmp(args_info.error_estimate_output_arg, "-"))
 176       ferr = stdout;
 177     else
 178       ferr = fopen(args_info.error_estimate_output_arg, "w");
 179     QPMS_ENSURE(ferr, "Could not open error output file %s",
 180         args_info.error_estimate_output_arg);
 181   if (args_info.output_given && !strcmp(args_info.output_arg, "-")
 182       && args_info.output_arg[0]) {
 183     out = fopen(args_info.output_arg, "w");
 184     QPMS_ENSURE(out, "Could not open output file %s", args_info.output_arg);
 185   } else
 186     out = stdout;
 187
 188   {
 189     const double omega = scuffomega * SCUFF_OMEGAUNIT;
 190     const double EeV = omega * hbar / eV;
 191     const double k0_vac = omega / c0;
 192     const double k0_eff  = k0_vac * refindex;
 193     const double eta = 5.224/a; // FIXME quite arbitrary, but this one should work
 194
 195     // indices : destpart (A/B-particle), srcpart (A/B-particle), coeff type (A/B- type), desty, srcy
 196     complex double W[npart][npart][2][c->nelem][c->nelem];
 197     double Werr[npart][npart][npart][c->nelem][c->nelem];
 198
 199     for (size_t ki = 0; ki < kcount; ++ki) {
 200       cart2_t beta = klist[ki];
 201       memset(W, 0, sizeof(W));
 202       if(ferr)
 203         memset(Werr, 0, sizeof(Werr));
 204
 205       const ptrdiff_t deststride = &(W[0][0][0][1][0]) - &(W[0][0][0][0][0]);
 206       const ptrdiff_t srcstride = &(W[0][0][0][0][1]) - &(W[0][0][0][0][0]);
 207       assert (srcstride == 1 && deststride == c->nelem);
 208
 209       for (size_t ps = 0; ps < npart; ++ps) for (size_t pd = 0; pd < npart; ++pd)
 210         // TODO optimize (calculate only once for each particle shift; especially if pd == ps)
 211         qpms_trans_calculator_get_AB_arrays_e32(c,
 212             &(W[pd][ps][0][0][0]), ferr ? &(Werr[pd][ps][0][0][0]) : NULL, // Adest, Aerr,
 213             &(W[pd][ps][1][0][0]), ferr ? &(Werr[pd][ps][1][0][0]) : NULL, // Bdest, Berr,
 214             deststride, srcstride,
 215             eta, k0_eff, b1, b2,
 216             beta,
 217             cart2_substract(part_positions[pd], part_positions[ps]),  // CHECKSIGN
 218             maxR, maxK
 219             );
 220         // TODO CHECK B<-A vs. A<-B relation
 221
 222       fprintf(out, "%.16g\t%.16g\t%.16g\t%.16g\t%.16g\t",
 223           scuffomega, EeV, k0_eff, beta.x, beta.y);
 224       if(ferr) fprintf(ferr, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
 225           scuffomega, EeV, k0_eff, beta.x, beta.y);
 226       size_t totalelems = sizeof(W) / sizeof(complex double);
 227       for (size_t i = 0; i < totalelems; ++i) {
 228         complex double w = ((complex double *)W)[i];
 229         fprintf(out, "%.16g\t%.16g\t", creal(w), cimag(w));
 230         if (ferr)
 231           fprintf(ferr, "%.3g\t", ((double *)Werr)[i]);
 232       }
 233       fputc('\n', out);
 234       if(ferr) fputc('\n', ferr);
 235     }
 236   }
 237
 238   fclose(out);
 239   if(ferr) fclose(ferr);
 240
 241 //#ifdef KSTDIN
 242   free(klist);
 243 //#endif
 244
 245   qpms_trans_calculator_free(c);
 246
 247 }