Merge branch 'trapping_debug' into 'master'

Trapping debug

See merge request giacomo.mulas/np_tmcode!13

build/README.md

@@ -62,6 +62,16 @@ The default behaviour of `np_sphere` is to take the same input files as `edfb` a
 
 ### trapping
 
-The execution of trapping programs requires at least one of the previous programs to have produced a complete output set. The working assumption is that the `frfme` program is a pre-requisite to `lffft`.
+The execution of trapping programs requires at least one of the previous programs to have produced a complete output set. A light-weight trapping calculation has been configured with input and legacy output files stored in the `../../test_data/trapping/` folder. Since the FORTRAN code assumes the input and output to be defined within the program, it is not yet possible to run the *FORTRAN* version on this case, unless the source code is modified accordingly. Conversely, the *C++* version can be executed without the need to modfy and re-compile the code. The work-flow to test trapping is described below.
 
-*TODO:* Define a common format for binary I/O operations on the TTMS file.
+1. cd to the `build/sphere` folder.
+2. run `np_sphere` with arguments to take input from the trapping tast and write output in the trapping folder:
+
+   > ./np_sphere ../../test_data/trapping/DEDFB ../../test_data/trapping/DSPH ../trapping
+
+3. cd to the trapping folder.
+4. run `np_trapping`
+
+   > ./np_trapping ../../test_data/trapping/DFRFME ../../test_data/trapping/DLFFT .
+
+5. Check the consistency between `np_trapping` output files (`c_out66.txt` and `c_out67.txt`) and the legacy *FORTRAN* output for this case (the files named, respectively, `fort.66` and `fort.67` in the `test_data/trapping/` folder).

src/libnptm/tra_subs.cpp

@@ -190,7 +190,7 @@ void ffrt(
 
   ctqce = new complex<double>[3]();
   ctqcs = new complex<double>[3]();
-  for (int l60 = 1; l60 < cil->le; l60++) {
+  for (int l60 = 1; l60 <= cil->le; l60++) {
     int lpo = l60 + 1;
     int il = l60 * lpo;
     int ltpo = l60 + lpo;
@@ -438,7 +438,6 @@ void wamff(
     un[0] = -sph;
     un[1] = cph;
     un[2] = 0.0;
-    // Would call PWMALP(W,UP,UN,YLM,LM)
     pwmalp(w, up, un, ylm, lm);
     double apfa = sth * apfafa;
     if (spd > 0.0) {
@@ -446,7 +445,7 @@ void wamff(
       double cthl = sqrt(1.0 - sthl * sthl);
       double arg = spd * (z - (r / rir) * cthl);
       cfam = (tra * std::exp(-apfa * apfa) / sqrt(cthl)) * std::exp(uim * arg);
-      if (lmode == 0) { // CHECK: Computed GO TO, not sure what happens with LMODE = 0
+      if (lmode == 0) {
 	if (sth == 0.0) { // label 45
 	  ftc1 = ftcn;
 	  ftc2 = ftcn;

src/sphere/sphere.cpp

@@ -294,10 +294,17 @@ void sphere(string config_file, string data_file, string output_path) {
 	    ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
 	    ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
 	    for (int lmi = 0; lmi < gconf->l_max; lmi++) {
-	      complex<double> element1 = -1.0 / c1->rmi[0][lmi];
-	      complex<double> element2 = -1.0 / c1->rei[0][lmi];
-	      ttms.write(reinterpret_cast<char *>(&element1), sizeof(complex<double>));
-	      ttms.write(reinterpret_cast<char *>(&element2), sizeof(complex<double>));
+	      complex<double> element1 = -1.0 / c1->rmi[lmi][0];
+	      complex<double> element2 = -1.0 / c1->rei[lmi][0];
+	      double vreal, vimag;
+	      vreal = element1.real();
+	      vimag = element1.imag();
+	      ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+	      ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
+	      vreal = element2.real();
+	      vimag = element2.imag();
+	      ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+	      ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
 	    }
 	    ttms.write(reinterpret_cast<char *>(&(sconf->radii_of_spheres[0])), sizeof(double));
 	    ttms.close();

src/trapping/clffft.cpp

@@ -217,7 +217,6 @@ void lffft(string data_file, string output_path) {
 	  // label 155
 	  if (!goto160) {
 	    if (jss != 1) {
-	      // Would open the ITT file.
 	      string tlfft_name = output_path + "/c_TLFFT";
 	      tlfft.open(tlfft_name.c_str(), ios::out | ios::binary);
 	      if (tlfft.is_open()) {
@@ -275,6 +274,7 @@ void lffft(string data_file, string output_path) {
 		      ws[ie] = wsl[iel];
 		    } // i175 loop
 		  }
+		} // i loop
 		// label 180
 		if (is != 2222) {
 		  if (is != 1111) {
@@ -371,7 +371,6 @@ void lffft(string data_file, string output_path) {
 		  delete[] ffte;
 		  delete[] ffts;
 		}
-		} // i loop
 	      } // ix475 loop
 	    } // iy475 loop
 	  } // iz475 loop

test_data/trapping/DEDFB

+    1  0
+ 1.765736D00  3.0D8   1.000000D00  0     1  0  3
+ 785.0000000D-09
+    1
+  1   0.200000D-06
+ 1.000000D00
+ (2.89D00,   0.03D00)
+0

test_data/trapping/DFRFME

-1  160
-0  8  25  10  10  10  10  1.0d0
-  1.2d0  2.0d0  1.0d0
-2.0d-05  0.15d0  2.3104d0
+1  16
+0  2  25  10  10  10  10  3.0d0
+  1.3d0  2.0d0  1.0d0
+0.0d-06  0.1733d0  2.3104d0
 e  1
 0
\ No newline at end of file

test_data/trapping/DSPH

+    1   2    0  149  300    0
+ 0.00D01  0.00D01  8.90D01  0.00D00  0.00D01  1.00D01
+ 0.00D01  0.00D01  4.90D01  0.00D00  0.00D00  1.00D01
+  1
+0
\ No newline at end of file

test_data/trapping/OEDFB

+  JXI     WLS          WNS          PUS          EVS          XIV
+    1   7.8500D-07   8.0041D+06   2.4012D+15   1.5805D+00   8.0041D+06
+  DIELECTRIC FUNCTIONS
+ SPHERE N.   1
+ NONTRANSITION LAYER N. 1, SCALE = WLS
+    1   2.8900D+00  3.0000D-02

test_data/trapping/OFRF

-  FSHMX = -7.705737229958699E-006
-  FRSH = -3.000000000000000E-006

test_data/trapping/OFRFME

+ IF JLML < NLMMT, PRESERVE TEMPTAPE1, TEMPTAPE2, AND TFRFME,
+ AND RESTART LM RUN WITH JLMF = JLML+1
+  FRSH =   0.0000000000000000     

test_data/trapping/OSPH

+ READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM
+     1    2    0  149  300    0
+ READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST
+  0.000D+00 0.000D+00 8.900D+01 0.000D+00 0.000D+00 1.000D+01
+ READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
+  0.000D+00 0.000D+00 4.900D+01 0.000D+00 0.000D+00 1.000D+01
+ READ(IR,*)JWTM
+     1
+  READ(ITIN)NSPHT
+  READ(ITIN)(IOG(I),I=1,NSPH)
+  READ(ITIN)EXDC,WP,XIP,IDFC,NXI
+  READ(ITIN)(XIV(I),I=1,NXI)
+  READ(ITIN)NSHL(I),ROS(I)
+  READ(ITIN)(RCF(I,NS),NS=1,NSH)
+ 
+  REFR. INDEX OF EXTERNAL MEDIUM=  1.3288100D+00
+   LIN
+ 
+========== JXI =  1 ====================
+  VK=  8.0040577D+06, XI=  8.0040577D+06
+     SPHERE  1
+  SIZE=  1.6008115D+00, REFRACTIVE INDEX=  1.7000229D+00  8.8234106D-03
+ ----- SCS ----- ABS ----- EXS ----- ALBEDS --
+  7.2050597D-14  5.9852332D-15  7.8035830D-14  9.2330147D-01
+ ---- SCS/GS -- ABS/GS -- EXS/GS ---
+  5.7336043D-01  4.7628972D-02  6.2098941D-01
+  FSAS=  1.6065129D-20  6.4776389D-21
+  QSCHU=  1.0965033D+00, PSCHU=  2.7194273D+00, S0MAG=  1.8784912D+00
+  COSAV=  6.7469268D-01, RAPRS=  2.9423820D-14
+  IPO= 1, TQEk=  0.0000000D+00, TQSk=  0.0000000D+00
+  IPO= 2, TQEk=  0.0000000D+00, TQSk=  0.0000000D+00
+********** JTH =  1, JPH =  1, JTHS =  1, JPHS =  1 ********************
+  TIDG= 0.000D+00, PIDG= 0.000D+00, TSDG= 0.000D+00, PSDG= 0.000D+00
+  SCAND= 0.000D+00
+  CFMP=  1.0000000D+00, SFMP=  0.0000000D+00
+  CFSP=  1.0000000D+00, SFSP=  0.0000000D+00
+  UNI=(-0.00000D+00, 1.00000D+00, 0.00000D+00)
+  UNS=(-0.00000D+00, 1.00000D+00, 0.00000D+00)
+     SPHERE  1
+  SAS(1,1)=  1.6065129D-20  6.4776389D-21, SAS(2,1)=  0.0000000D+00 -0.0000000D+00
+  SAS(1,2)=  0.0000000D+00 -0.0000000D+00, SAS(2,2)=  1.6065129D-20  6.4776389D-21
+  Fx=  0.0000000D+00, Fy=  0.0000000D+00, Fz=  2.9423820D-14
+  MULS
+          3.1194076D-14  0.0000000D+00 -0.0000000D+00 -0.0000000D+00
+          0.0000000D+00  3.1194076D-14  0.0000000D+00  0.0000000D+00
+         -0.0000000D+00  0.0000000D+00  3.1194076D-14  0.0000000D+00
+         -0.0000000D+00  0.0000000D+00  0.0000000D+00  3.1194076D-14
+  MULSLR
+          3.1194076D-14  0.0000000D+00 -0.0000000D+00 -0.0000000D+00
+          0.0000000D+00  3.1194076D-14 -0.0000000D+00 -0.0000000D+00
+         -0.0000000D+00 -0.0000000D+00  3.1194076D-14  0.0000000D+00
+         -0.0000000D+00 -0.0000000D+00  0.0000000D+00  3.1194076D-14