Remove vacancy limit for EADL relaxations

Remove statements in routine egs_init_relax to allow all shells with
binding energy above 1 keV available in the Evaluated Atomic Data
Library (EADL). Also remove some commented out portions of the code.

Up to now, the default relaxation implementation in EGSnrc allowed
vacancies from atomic shells K and L1-L3, as well as <M> and <N> shells
in an averaged way. Also, different interactions consider different sets
of shells:

EGSnrc shell types:
-------------------
Compton interaction (shn_array)  1..7: K, L1, L2, L3, <M>, <N>, rest
Photoelectric effect             1..5: K, L1, L2, L3, <M>
Electron impact ionization       1..4: K, L1, L2, L3

However, the EADL relaxation implementation can actually accept more
shell types, up to what is available in the library (27 shells for
atomic number Z=99). But EGSnrc was forcing M1 and N1 shells instead of
<M> and <N>, and using O1 as the outermost shell, via the following
statements in egs_init_relax:

IF (shell_num(ish)<=7) shell_eadl(iZ,shell_num(ish)) = ish;
...
shell_eadl(iZ,6) = shell_ntot + 10; "reseting to type 16 (N1)"
shell_eadl(iZ,7) = shell_ntot + 17; "reseting to type 27 (O1)"

Since EGSnrc shell types can range from 1 to 7, vacancies can
potentially be created in M1, N1, or O1 shells if energetically
possible.

PROBLEM: This approach ignores that M2 and M3 vacancies can in fact be
favored energetically! For instance for Z=74, M1 and M2 vacancies with
binding energies 2.8 keV and 2.6 keV respectively could generate
characteristic x-ray lines above 1 keV.

SOLUTION: Remove the aforementioned statements in egs_init_relax,
allowing all shells with binding energy above 1 keV.

OUTCOME: Invoking relax using shells 5, 6, or 7 will create vacancies in
the M1, M2, and M3 shells. This is not entirely correct for the
photoelectric effect, where the sampling is done for <M> only, nor for
Compton since shell number 5 will then correspond to shells M1 to M5,
and shell number 6 to shells N1 to N7. For Compton this can be corrected
by using shell_array instead of shn_array in the call to relax.

HEN_HOUSE/src/egsnrc.mortran

@@ -5471,15 +5471,9 @@ DO medio = 1,nmed [
         curr_rec  = curr_rec+1;
         read(relax_unit, rec=curr_rec) be_r; shell_be(ish) = be_r;
 
-        /*
-           IF (iZ=74)[
-           write(*,*) 'Shell# ', ish ,' of type ', shell_type(ish);
-           write(*,*) ' -> be = ', shell_be(ish) ,' MeV ';
-           write(*,*) 'has ', ntran ,' transitions ';
-       ] */
         shell_Z(ish) = iZ;
         shell_num(ish) = ish - shell_ntot;
-        IF (shell_num(ish)<=7) shell_eadl(iZ,shell_num(ish)) = ish;
+        shell_eadl(iZ,shell_num(ish)) = ish;
         DO k=1,ntran[
             curr_rec  = curr_rec+1;
             read(relax_unit,rec=curr_rec) itmp(k);
@@ -5490,10 +5484,6 @@ DO medio = 1,nmed [
               radiative and by 65 for non-radiative transitions*/
             IF (itmp(k)<64) [itmp(k) +=1;]
             ELSE            [itmp(k) += 65;]
-          /*
-             IF (iZ=74)[
-             write(*,*) ' transition ',k,': type=',itmp(k),' prob=',wtmp(k);
-          ]*/
         ]
         IF( shell_be(ish) < min_be ) [
             relax_first(ish) = -1;
@@ -5523,26 +5513,6 @@ DO medio = 1,nmed [
             relax_ntot = relax_ntot + ntran;
         ]
       ]
-      shell_eadl(iZ,6) = shell_ntot + 10; "reseting to type 10 (N1)"
-      shell_eadl(iZ,7) = shell_ntot + 17; "reseting to type 27 (O1)"
-      /*
-      IF (iZ=74)[
-       $egs_info('(a,i4)','For Z = 74: ish_ini = ',shell_ntot+1);
-          sh_eadl = shell_eadl(iZ,1);
-          $egs_info('(4(a,i3))','shell_egs = ',1,
-                  ' shell_eadl = ',sh_eadl,
-                  ' shell type = ',shell_type(sh_eadl),
-                  ' # of transitions = ', relax_ntran(shell_ntot+1));
-          ntran = relax_ntran(shell_ntot+1);
-          IF (sh_eadl>0)[
-            DO m=0,ntran-1[
-              $egs_info(*,' transition ', m+1,': ',
-                     relax_state(relax_first(shell_ntot+1) + m),',',
-                     relax_prob( relax_first(shell_ntot+1)+ m));
-            ]
-          ]
-       $egs_info('(a,i4)','ish_final = ',shell_ntot+nshell);
-      ]*/
 
       shell_ntot = shell_ntot + nshell;
 
@@ -5604,16 +5574,10 @@ Nvac_max = -1;
 shell = shell_eadl(iZ,shell_egs);
 IF( shell < 1 | shell > $MAXSHELL ) [ return; ] "unknown vacancy"
 
-"$egs_info('(3(a,i3))','Z = ',iZ,"
-"        ' shell_egs = ',shell_egs,' shell_eadl = ', shell);"
-
 irl = ir(np);
 Ec = ecut(irl) - rm; Pc = pcut(irl);
 min_E = min(Ec,Pc); min_E = max($RELAX-CUTOFF,min_E);
 Evac = shell_be(shell);
-"IF( Evac < min_E | relax_ntran(shell) < 1 ) ["
-"    edep += Evac; return;"
-"]"
 rfu_Z = shell_Z(shell); rfu_j0 = shell;
 rfu_n0 = shell_num(shell); rfu_t0 = shell_type(shell);
 rfu_E0 = Evac;
@@ -5638,35 +5602,22 @@ LOOP [
         new_state = relax_state(relax_first(vac)+new_state-1);
      IF( new_state <= 64 ) [ "It was a radiative transition"
         iqf = 0; new_state += vac - shell_num(vac); Ef = shell_be(new_state);
-        "IF( Ef > min_E & relax_ntran(new_state) > 0 ) [ "
         Nvac += 1; vacs(Nvac) = new_state;
-            "IF (Nvac>Nvac_max) Nvac_max = Nvac;"
-        "]"
-        "ELSE [ edep += Ef; NEXT;]"
         Ecc = Pc;
-"IF (Ef>Evac)[$egs_info('(2(a,i4))',"
-"'vac_type=',shell_type(vac),' trans_type=',shell_type(new_state));]"
      ]
      ELSE [ "It was a non-radiative transition (Auger or Coster-Kronig)"
         iqf = -1; new1 = new_state/64; new2 = new_state - 64*new1;
         new1 += vac - shell_num(vac); new2 += vac - shell_num(vac);
         Ef1 = shell_be(new1); Ef2 = shell_be(new2);
-        "IF( Ef1 > min_E & relax_ntran(new1) > 0 ) [ "
         Nvac += 1; vacs(Nvac) = new1;
-            "IF (Nvac>Nvac_max) Nvac_max = Nvac;"
-        "] ELSE [ edep += Ef1;]"
-        "IF( Ef2 > min_E & relax_ntran(new2) > 0 ) [ "
         Nvac += 1; vacs(Nvac) = new2;
-            "F (Nvac>Nvac_max) Nvac_max = Nvac; "
-        "] ELSE [ edep += Ef2; NEXT;]"
         Ef = Ef1 + Ef2; Ecc = Ec;
      ]
     ]
     Ex = Evac - Ef;
     edep_local = 0;
     IF( Ex > Ecc ) [
         np += 1; $CHECK-RELAX-STACK(np,'new_relax');
-        "$egs_info('(a,i9,a,f7.3)','-> np=',np,' iq=',iqf);"
         iq(np) = iqf;
         $TRANSFER PROPERTIES TO (np) FROM (np_save);
         $RANDOMSET rnno; cost = 2*rnno-1; sint = 1-cost*cost;
@@ -5683,7 +5634,6 @@ LOOP [
         ELSE [
             e(np) = Ex + rm; $AUSCALL($AUGERTRA);
         ]
-        "$egs_info('(a,i2,a,f7.3,a)','-> iq=',iqf,' E=',Ex*1000.,' keV');"
     ] ELSE [
         edep += Ex;
 
@@ -5699,7 +5649,6 @@ LOOP [
     IF( Nvac = 0 ) EXIT;
     vac = vacs(Nvac); Evac = shell_be(vac); Nvac -= 1;
 ]
-"IF (Nvac_max>1) $egs_info('(a,i2,a)','-> ',Nvac_max,' vacancies');"
 return;
 end;
 ;