Fix #749: add SAD to z position for phsp scoring

Add the source-axis distance to the z position of the particles before
the phase space is scored. This change only affects IAEA phsp sources,
and only when scoring particles in BEAMnrc coordinates.

HEN_HOUSE/user_codes/dosxyznrc/dosxyznrc.mortran

@@ -3727,6 +3727,11 @@ IF(iarg=5 & i_phsp_out~=0 & ir(np)=irmax+1 & irold<ir(np))[
  w(np)=r_13(k_field)*u_temp + r_23(k_field)*v_temp +
  r_33(k_field)*w_temp;
 
+ IF(i_iaea_in=1)[
+ "for IAEA phsp source we need to add SAD back to Z coordinate"
+ z(np) = z(np) + iaea_SAD;
+ ]
+
  "Write the particle to a phsp"
  $DOSXYZ_WRITE_PHSP;
 

HEN_HOUSE/user_codes/dosxyznrc/srcxyznrc.macros

@@ -166,6 +166,7 @@ REPLACE {;COMIN/SOURCE/;} WITH {
  esrc_sp, "kinetic energy of source particle (enflag=1) "
  SSD, "Source to surface dist for point source particles ( isource=3)"
  einsrcold,"stores einsrc for source 20 (isource=20) "
+ iaea_SAD, "stores SAD for IAEA phase space sources (isource=2,8,20)"
  NINCSRC, "No. of particles incident from original source (isource=2,8) "
  r_dbs, "DBS splitting radius used to generate phsp source (isource=2,8)"
  ssd_dbs, "SSD at which above DBS radius defined (isource=2,8)"
@@ -243,7 +244,7 @@ REPLACE {;COMIN/SOURCE/;} WITH {
  pang,angfixed,angmin,angmax,pgang,r_11, r_12, r_13, r_21,
  r_22, r_23, r_31, r_32, r_33,dsource,dsourcetemp,
  ein, eksrcm, esrc_sp, SSD, NINCSRC,
- r_dbs, ssd_dbs,z_dbs,survival_ratio,einsrcold;
+ r_dbs, ssd_dbs,z_dbs,survival_ratio,einsrcold,iaea_SAD;
  $INTEGER isource, iqin, iqphsp, iqinc,latchold,iqinold,MLCflag,SHLflag,
  calflag, ixinl, ixinu, jyinl, jyinu, kzinl, kzinu,
  klowx, khix, klowy, khiy, klowz, khiz,

HEN_HOUSE/user_codes/dosxyznrc/srcxyznrc.mortran

@@ -3384,7 +3384,6 @@ $REAL xsrc,ysrc,zsrc, "X,Y,Z position of incident particle (isource=2,8,9)"
  uin,vin,win, "U,V,W of incident particle"
  enin, "energy of incident particle"
  zsrcdum, "dummy value of zsrc read in from BEAM sim (isource=9)"
- SAD_s20, "for source 20 if iaea phsp is read SAD is needed"
  r1,rnno1,rnno2, "random no.s"
  fw,
  rxyz, "radial position of incident particle"
@@ -3408,6 +3407,10 @@ $LONG_INT n_hist_dum;
 save wsrc; "in case of phsp particle recycling"
 
 zsrc = dsource;
+iaea_SAD = 0;
+IF(i_iaea_in = 1)[
+ iaea_SAD = -dsource; "save SAD for iaea phsp sources"
+]
 
 "----------------------------------------------------------------"
 IF(isource = 0)[ "Frontal parallel beam source"
@@ -4153,8 +4156,7 @@ IF(NofREPEAT > 0 & ISMOOTH = 1)[
 
 "Rotate/translate the source particle into position
 IF(isource<7 | isource=9)[
- IF(i_iaea_in=1) zsrc = zsrc + dsource;
- "this is zsrc - SAD, but dsource has been set -ve in srcinit"
+ IF(i_iaea_in=1) zsrc = zsrc - iaea_SAD;
 
  uin = r_11(1)*usrc + r_12(1)*vsrc + r_13(1)*wsrc;
  vin = r_21(1)*usrc + r_22(1)*vsrc + r_23(1)*wsrc;
@@ -4180,9 +4182,9 @@ ELSEIF(isource=20 | isource=21)[
  IF(i_iaea_in=1)[
  "for iaea phase space z is read from the particle"
  " dsourcetemp is now supposed to be SAD"
- SAD_s20 = (dsourcetemp(K-1)+(dsourcetemp(K)-dsourcetemp(K-1))*
+ iaea_SAD = (dsourcetemp(K-1)+(dsourcetemp(K)-dsourcetemp(K-1))*
  (frMU_indx-muIndex(K-1))/(muIndex(K)-muIndex(K-1)));
- zsrc = zsrc - SAD_s20;
+ zsrc = zsrc - iaea_SAD;
  ] ELSE [
  zsrc = -1.0*(dsourcetemp(K-1)+(dsourcetemp(K)-dsourcetemp(K-1))*
  (frMU_indx-muIndex(K-1))/(muIndex(K)-muIndex(K-1)));
@@ -4229,6 +4231,7 @@ cosphi = cos(phi(1)*3.141593/180.); sinphi = sin(phi(1)*3.141593/180.);
 
 ]
 ELSE[ "choose an incident angle from the distribution"
+ IF(i_iaea_in=1) zsrc = zsrc - iaea_SAD; "in case we have IAEA isource=8"
  $RANDOMSET RNNO1;
  k_field=numang;"use global variable k_field so we can untransform"
  "if scoring phsp in BEAM coordinate system"