# extract_flux.r
# 3D flux files should be present under ./flux/flux_Arplus_*.txt

integer IZ0 = 13,   IZ1 = 86;   # Mininum and maximum Z coordinate index
float   XM  = 167.5, YM = 247;  # Center coordinate of cylinder in XY plane
string BASENAME : "flux_Arplus_$(time%,2f)us";  # Scheme of input filenames  

integer ColX=1, ColY=2, ColZ=0; # Column offsets for rotating the coordinate system
                                # If the cylindricat target in the geo file is already
                                # oriented in Z direction, use ColX=0, ColY=1, ColZ=2

# Number of entries in X, Y, Z
# Lowest and highest coordinate in X, Y, Z, respectively

integer NX = 70;  float X0 = 47.03214285714284, X1 = 258.4678571428572, DX = (X1-X0)/(NX-1);
integer NY = 37;  float Y0 = 267,               Y1 = 339,               DY = (Y1-Y0)/(NY-1);
integer NZ = 100; float Z0 = 3.75,              Z1 = 746.25,            DZ = (Z1-Z0)/(NZ-1);
 
float Alpha0 = -30, Alpha1=50, dAlpha=1;  # Extracted angular range
float R0     = 69;                        # Target diameter plus safety margin
float DR     = 2;                         # Safety margin in the order of at least one half cell spacing

# =======================================================================================================
# No more adjustments needed from this point on
# =======================================================================================================

string str12(x=float) = return "$("$(x)"%12s)";
integer dim_fn = (IZ1-IZ0+1)*((Alpha1-Alpha0)/dAlpha+1);
float FN_Sum[dim_fn];

float realX(x=float, y=float, z=float) = return ColX==0 ? x ! ColY==0 ? y ! z;
float realY(x=float, y=float, z=float) = return ColX==1 ? x ! ColY==1 ? y ! z;
float realZ(x=float, y=float, z=float) = return ColX==2 ? x ! ColY==2 ? y ! z;

block extract_flux(BASENAME=string, sum->var) : {
  float bilinear(f00=float, f10=float, f01=float, f11=float, x=float, y=float, x0=float, y0=float, x1=float, y1=float) = {
    float F = 1/((x1-x0)*(y1-y0));
    return F*(f00*(x1-x)*(y1-y) + f10*(x-x0)*(y1-y) + 
	      f01*(x1-x)*(y-y0) + f11*(x-x0)*(y-y0));
  };

  float index(x=float, y=float, z=float) = 
    return floor((x-X0)/DX) + NX*floor((y-Y0)/DY) + NX*NY*floor((z-Z0)/DZ);

  float simeq(a=float, b=float) = return abs(a-b)<1e-6 ? 1 ! 0;

  console "** Reading $(BASENAME).txt **\n";
  SHEET M={
    filename="flux/$(BASENAME).txt";
    startline=1; separator=" ";
  };


  integer iz;
  float   zplane -> Z0 + DZ*iz;

  console "** Generating flux fields **\n";
  float FX[NX*NY*NZ], FY[NX*NY*NZ], FZ[NX*NY*NZ];

  integer i, ind;
  for(i=0; i<M.n_rows; i=i+1) {
    float x = floor((M.data[4+ColX][i]-X0)/DX+0.5)*DX+X0+DX/10,
          y = floor((M.data[4+ColY][i]-Y0)/DY+0.5)*DY+Y0+DY/10,
          z = floor((M.data[4+ColZ][i]-Z0)/DZ+0.5)*DZ+Z0+DZ/10;

    float fx= M.data[7+ColX][i],
	  fy= M.data[7+ColY][i],
	  fz= M.data[7+ColZ][i];

    ind = index(x, y, z);  
    FX[ind] = fx; FY[ind] = fy; FZ[ind] = fz;
  }  

  for(iz=IZ0; iz<=IZ1; iz=iz+1) {
    console "** Processing plane ix=$(iz) **\n";

    if(iz==IZ0) then file "$(BASENAME)_EXT.pos" << out "View\"$(BASENAME)\" {\n";

    append "$(BASENAME)_EXT.pos" << {
      float alpha, pi = 4*atan(1), x, y, x0, y0, z0, nx, ny, tx, ty;
      integer ind;

      for(alpha=Alpha0; alpha<=Alpha1+1e-6; alpha=alpha+dAlpha) {
	x =  XM+(R0+DR)*sin(alpha*pi/180);
	y =  YM+(R0+DR)*cos(alpha*pi/180);

	x0 = floor((x-X0)/DX)*DX + X0;
	y0 = floor((y-Y0)/DY)*DY + Y0;

	ind = index(x, y, zplane+1e-3);

        fx = bilinear(FX[ind], FX[ind+1], FX[ind+NX], FX[ind+NX+1], x, y, x0, y0, x0+DX, y0+DY);
        fy = bilinear(FY[ind], FY[ind+1], FY[ind+NX], FY[ind+NX+1], x, y, x0, y0, x0+DX, y0+DY);
        fz = bilinear(FZ[ind], FZ[ind+1], FZ[ind+NX], FZ[ind+NX+1], x, y, x0, y0, x0+DX, y0+DY);

	# Normal unit vector
        nx = (x-XM)/sqrt((x-XM)^2+(y-YM)^2);
        ny = (y-YM)/sqrt((x-XM)^2+(y-YM)^2);

	# Tangential unit vector
	tx = -ny;
	ty = nx;

	fn = fx*nx + fy*ny;    # Normal flux component
	ft = fx*tx + fy*ty;    # Tangential flux component

	out "SP($(realX(x, y, zplane)), $(realY(x, y, zplane)), $(realZ(x, y, zplane))){$(-fn)};\n";

        # fn auf externes Array addieren
        integer isum = (alpha-Alpha0)/dAlpha + ((Alpha1-Alpha0)/dAlpha+1)*(iz-IZ0);
        sum[isum] = sum[isum] + fn;
      }
    };

    if(iz==IZ1) then append "$(BASENAME)_EXT.pos" << out "};\n";
  }


  return 0;
};

integer count;
for(time=50.0; time<=60.0+1e-6; time=time+1) {
  extract_flux(BASENAME, FN_Sum);
  count=count+1;
}

for(i=0; i<FN_Sum.len; i=i+1) FN_Sum[i] = -(FN_Sum[i])/count;  # Flux unit is 1/m^2s

# Output of averaged flux as POS and TXT
file "AVG_flux_Arplus.pos" << {
  file "AVG_flux_Arplus.txt" << out "$("alpha"%12s)  $("z"%12s)  $("fn"%12s)\n";
  out "View\"AVG_Flux_Arplus\" {\n";
  integer iz;
  float   zplane -> Z0 + DZ*iz;
  for(iz=IZ0; iz<=IZ1; iz=iz+1) {
    float alpha, pi = 4*atan(1), x, y, x0, y0, z0, nx, ny, tx, ty;
    for(alpha=Alpha0; alpha<=Alpha1+1e-6; alpha=alpha+dAlpha) {
      x =  XM+(R0+DR)*sin(alpha*pi/180);
      y =  YM+(R0+DR)*cos(alpha*pi/180);
      integer isum = (alpha-Alpha0)/dAlpha + ((Alpha1-Alpha0)/dAlpha+1)*(iz-IZ0);
      out "SP($(realX(x, y, zplane)), $(realY(x, y, zplane)), $(realZ(x, y, zplane))){$FN_Sum[isum]};\n";
      append "AVG_flux_Arplus.txt" << out "$(str12(alpha))  $(str12(zplane))  $(str12(FN_Sum[isum]))\n";
    }
  }
  out "};\n";
};