/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meq.h"
# ifdef SINGLE
# define FCT  svizr
# else
# define FCT  dvizr
# endif
/* function [Wpol,Ft0,Vp] = vizr(Fx,lyp,ry,iry,rBt,drx,dzx,nrx,nzx) */
void FCT(FLT *y1, FLT *y2, FLT *y3, FLT *f, bool *l, FLT*r, FLT *ir, FLT rBt, 
								 FLT dr, FLT dz, int nr, int nz) {
	int i,j,ni;
	FLT x, yr, yz, yp = 0.0, yi = 0.0, yv = 0.0,
     drdz = dr/dz, dzdr = dz/dr, ds = dr*dz,
     *f1 = f+1, *f2 = f+nz, *f3 = f2+2, *f4 = f2+(nz+1);
 for (j=nr-2;j--;) {
  yr = yz = 0.0;
  ni = 0; // counter of how many grid points on this R contain plasma
  for (i=nz-2;i--;) {
   x = *f4++ - *f1++;
   yr += *l * x * x;
   x = *f3++ - *f2++;
   yz += *l * x * x;
   ni += *l++;
  }
  yp += (yr*dzdr + yz*drdz) * *ir;
  yi += ni * *ir++;
  yv += ni * (*r++);
  f1 += 2; f2 += 2; f3 += 2; f4 += 2;
 }
 *y1 = 15831.43494411528*yp; /* 1/(16*pi*mu0) */
 *y2 = rBt*ds*yi;
 *y3 = 6.283185307179586*ds*yv;  // Volume = 2pi\int RdRdZ = 2*pi*sum(ni*Ri)*dr*dz
 return;
}
