diff --git a/.gitignore b/.gitignore
index 55ffb472..6c3d0012 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,10 @@
-# This is currently empty. If you need some ignores, which are specific to your workflow (e.g. IDE), consider
+# If you need some ignores, which are specific to your workflow (e.g. IDE), consider
 # including them into $GIT_DIR/info/exclude (in this repository) or globally into $XDG_CONFIG_HOME/git/ignore
 # See https://git-scm.com/docs/gitignore 
+
+# for MacOS
 .DS_Store
+
+# for Eclipse
+.сproject
+.project 
diff --git a/src/CalculateE.c b/src/CalculateE.c
index 2bd43163..c6c41a87 100644
--- a/src/CalculateE.c
+++ b/src/CalculateE.c
@@ -44,6 +44,7 @@ extern const Parms_1D phi_sg;
 extern const double ezLab[3],exSP[3];
 // defined and initialized in GenerateB.c
 extern const double C0dipole,C0dipole_refl;
+extern int vorticity;
 // defined and initialized in param.c
 extern const bool store_int_field,store_dip_pol,store_beam,store_scat_grid,calc_Cext,calc_Cabs,
 	calc_Csca,calc_vec,calc_asym,calc_mat_force,store_force,store_ampl;
@@ -172,6 +173,7 @@ void MuellerMatrix(void)
 	double matrix[4][4];
 	double theta,phi,ph,max_err,max_err_c2,max_err_s2,max_err_c4,max_err_s4;
 	doublecomplex s1,s2,s3,s4;
+	doublecomplex vort;
 	char fname[MAX_FNAME];
 	int i;
 	size_t index,index1,k_or,j,n,ind;
@@ -197,6 +199,11 @@ void MuellerMatrix(void)
 				si=sin(alph);
 				for (i=0;i<nTheta;i++) {
 					// transform amplitude matrix, multiplying by rotation matrix (-alpha)
+					/* Note, that amplitude matrix for vortex beams (e.g., Bessel ones) have to be additionally
+					 * multiplied by the phase factor exp(-I*vorticity*alpha) and additional factor I^vorticity for
+					 * scat_plane (see below). However, it does not change the Mueller matrix and we do not save the
+					 * amplitude matrix here. Hence, these phase factors are ignored for orientation averaging.
+					 */
 					if (yzplane) { // here the default (alpha=0) is yz-plane, so par=Y, per=X
 						s2 =  co*ampl_alphaY[index+1] + si*ampl_alphaX[index+1]; // s2 =  co*s20 + si*s30
 						s3 = -si*ampl_alphaY[index+1] + co*ampl_alphaX[index+1]; // s3 = -si*s20 + co*s30
@@ -243,15 +250,19 @@ void MuellerMatrix(void)
 				FCloseErr(mueller,F_MUEL,ONE_POS);
 			}
 		}
-		if (scat_plane) { // par=X, per=-Y
+		if (scat_plane) { // par=X*i^n, per=-Y*i^n (n - vorticity)
 			if (store_ampl) {
 				SnprintfErr(ONE_POS,fname,MAX_FNAME,"%s/"F_AMPL,directory);
 				ampl=FOpenErr(fname,"w",ONE_POS);
 				fprintf(ampl,THETA_HEADER" "AMPL_HEADER"\n");
+				vort=cpow(I,vorticity);
 				for (i=0;i<nTheta;i++) {
 					theta=i*dtheta_deg;
-					fprintf(ampl,ANGLE_FORMAT" "AMPL_FORMAT"\n",theta,REIM(-EplaneY[2*i]),REIM(EplaneX[2*i+1]),
-						REIM(-EplaneY[2*i+1]),REIM(EplaneX[2*i]));
+					s1=-EplaneY[2*i]*vort;
+					s2=EplaneX[2*i+1]*vort;
+					s3=-EplaneY[2*i+1]*vort;
+					s4=EplaneX[2*i]*vort;
+					fprintf(ampl,ANGLE_FORMAT" "AMPL_FORMAT"\n",theta,REIM(s1),REIM(s2),REIM(s3),REIM(s4));
 				}
 				FCloseErr(ampl,F_AMPL,ONE_POS);
 			}
@@ -261,6 +272,7 @@ void MuellerMatrix(void)
 				fprintf(mueller,THETA_HEADER" "MUEL_HEADER"\n");
 				for (i=0;i<nTheta;i++) {
 					theta=i*dtheta_deg;
+					// here we ignore common phase factor in the amplitude matrix (vort above)
 					ComputeMuellerMatrix(matrix,-EplaneY[2*i],EplaneX[2*i+1],-EplaneY[2*i+1],EplaneX[2*i],theta);
 					fprintf(mueller,ANGLE_FORMAT" "MUEL_FORMAT"\n",theta,COMP44M(matrix));
 				}
@@ -272,10 +284,11 @@ void MuellerMatrix(void)
 			/* compute Mueller Matrix in full space angle.
 			 * E-fields are stored in arrays EgridX and EgridY for incoming X and Y polarized light.
 			 * It is converted to the scattering matrix elements (see e.g Bohren and Huffman) :
-			 * s2 = cos(phi)E'X'par + sin(phi)E'Y'par
-			 * s3 = sin(phi)E'X'par - cos(phi)E'Y'par
-			 * s4 = cos(phi)E'X'per + sin(phi)E'Y'per
-			 * s1 = sin(phi)E'X'per - cos(phi)E'Y'per
+			 * s2 = vort*(cos(phi)E'X'par + sin(phi)E'Y'par)
+			 * s3 = vort*(sin(phi)E'X'par - cos(phi)E'Y'par)
+			 * s4 = vort*(cos(phi)E'X'per + sin(phi)E'Y'per)
+			 * s1 = vort*(sin(phi)E'X'per - cos(phi)E'Y'per)
+			 * where additional phase factor vort=exp(I*vorticity*(PI_OVER_TWO-ph)) is for vortex beams;
 			 * from these the mueller matrix elements are computed
 			 */
 			// open files for writing
@@ -330,15 +343,16 @@ void MuellerMatrix(void)
 					if (angles.type==SG_GRID) phi=angles.phi.val[j];
 					else phi=angles.phi.val[ind]; // angles.type==SG_PAIRS
 					ph=Deg2Rad(phi);
+					vort=imExp(vorticity*(PI_OVER_TWO-ph));
 					// rather complicated (but general) approach to determine rotation angle from XY to scattering plane
 					SetScatPlane(cthet,sthet,ph,tmp3,incPolper);
 					co=-DotProd(incPolper,incPolY);
 					si=DotProd(incPolper,incPolX);
 					// transform the amplitude matrix, multiplying by rotation matrix from per-par to X-Y
-					s2 = co*EgridX[index+1] + si*EgridY[index+1]; // s2 =  co*s20 + si*s30
-					s3 = si*EgridX[index+1] - co*EgridY[index+1]; // s3 =  si*s20 - co*s30
-					s4 = co*EgridX[index] + si*EgridY[index]; // s4 =  co*s40 + si*s10
-					s1 = si*EgridX[index] - co*EgridY[index]; // s1 =  si*s40 - co*s10
+					s2 = vort*(co*EgridX[index+1] + si*EgridY[index+1]); // s2 =  co*s20 + si*s30
+					s3 = vort*(si*EgridX[index+1] - co*EgridY[index+1]); // s3 =  si*s20 - co*s30
+					s4 = vort*(co*EgridX[index] + si*EgridY[index]); // s4 =  co*s40 + si*s10
+					s1 = vort*(si*EgridX[index] - co*EgridY[index]); // s1 =  si*s40 - co*s10
 					index+=2;
 
 					if (store_mueller) ComputeMuellerMatrix(matrix,s1,s2,s3,s4,theta);
diff --git a/src/GenerateB.c b/src/GenerateB.c
index 9a818851..577622bd 100644
--- a/src/GenerateB.c
+++ b/src/GenerateB.c
@@ -46,19 +46,7 @@ extern const opt_index opt_beam;
 
 // used in CalculateE.c
 double C0dipole,C0dipole_refl; // inherent cross sections of exciting dipole (in free space and addition due to surface)
-
-// used in crosssec.c
-double beam_center_0[3]; // position of the beam center in laboratory reference frame
-/* complex wave amplitudes of secondary waves (with phase relative to particle center);
- * The transmitted wave can be inhomogeneous wave (when msub is complex), then eIncTran (e) is normalized
- * counter-intuitively. Before multiplying by tc, it satisfies (e,e)=1!=||e||^2. This normalization is consistent with
- * used formulae for transmission coefficients. So this transmission coefficient is not (generally) equal to the ratio
- * of amplitudes of the electric fields. In particular, when E=E0*e, ||E||!=|E0|*||e||, where
- * ||e||^2=(e,e*)=|e_x|^2+|e_y|^2+|e_z|^2=1
- *
- * !!! TODO: determine whether they are actually needed in crosssec.c, or make them static here
- */
-doublecomplex eIncRefl[3],eIncTran[3];
+int vorticity;                 // Vorticity of vortex beams (besN for Bessel beams)
 // used in param.c
 const char *beam_descr; // string for log file with beam parameters
 
@@ -68,11 +56,23 @@ static double scale_x,scale_z; // multipliers for scaling coordinates
 static doublecomplex ki,kt;    // abs of normal components of k_inc/k0, and ktran/k0
 static doublecomplex ktVec[3]; // k_tran/k0
 static double p0;              // amplitude of the incident dipole moment
+#ifndef NO_FORTRAN
+static int besN;                  // Bessel beam order
+static double besAlpha;           // half-cone angle (in radians)
+static doublecomplex besKt,besKz; // wave-vector components (transverse and longitudinal)
+static doublecomplex besM[4];     // Matrix M defining the generalized Bessel beam
+#endif
 /* TO ADD NEW BEAM
  * Add here all internal variables (beam parameters), which you initialize in InitBeam() and use in GenerateB()
  * afterwards. If you need local, intermediate variables, put them into the beginning of the corresponding function.
  * Add descriptive comments, use 'static'.
  */
+ 
+// EXTERNAL FUNCTIONS
+
+#ifndef NO_FORTRAN 
+void bjndd_(const int *n,const double *x,double *bj,double *dj,double *fj);
+#endif
 
 //======================================================================================================================
 
@@ -84,15 +84,20 @@ void InitBeam(void)
 	/* TO ADD NEW BEAM
 	 * Add here all intermediate variables, which are used only inside this function.
 	 */
-
+	
 	// initialization of global option index for error messages
 	opt=opt_beam;
 	// beam initialization
+	beam_asym=(beam_center_0[0]!=0 || beam_center_0[1]!=0 || beam_center_0[2]!=0);
+	if (!beam_asym) vInit(beam_center); // not to calculate it for each orientation
+	vorticity = 0;
 	switch (beamtype) {
 		case B_PLANE:
 			if (IFROOT) beam_descr="plane wave";
-			beam_asym=false;
 			if (surface) {
+				/* here we assume that prop_0 will not change further (e.g., by rotation of particle), 
+				 * i.e. prop=prop_0 in GenerateBeam() below
+				 */				
 				if (prop_0[2]==0) PrintError("Ambiguous setting of beam propagating along the surface. Please specify "
 					"the incident direction to have (arbitrary) small positive or negative z-component");
 				if (msubInf && prop_0[2]>0) PrintError("Perfectly reflecting surface ('-surf ... inf') is incompatible "
@@ -136,21 +141,17 @@ void InitBeam(void)
 			}
 			return;
 		case B_DIPOLE:
-			vCopy(beam_pars,beam_center_0);
 			if (surface) {
 				if (beam_center_0[2]<=-hsub)
 					PrintErrorHelp("External dipole should be placed strictly above the surface");
 				inc_scale=1; // but scaling of Mueller matrix is weird anyway
 			}
-			// in weird scenarios the dipole can be positioned exactly at the origin; reused code from Gaussian beams
-			beam_asym=(beam_center_0[0]!=0 || beam_center_0[1]!=0 || beam_center_0[2]!=0);
-			if (!beam_asym) vInit(beam_center);
 			/* definition of p0 is important for scaling of many scattering quantities (that are normalized to incident
 			 * irradiance). Alternative definition is p0=1, but then the results will scale with unit of length
 			 * (breaking scale invariance)
 			 */
 			p0=1/(WaveNum*WaveNum*WaveNum);
-			if (IFROOT) beam_descr=dyn_sprintf("point dipole at "GFORMDEF3V,COMP3V(beam_center_0));
+			if (IFROOT) beam_descr="point dipole";
 			return;
 		case B_LMINUS:
 		case B_DAVIS3:
@@ -159,9 +160,6 @@ void InitBeam(void)
 			// initialize parameters
 			w0=beam_pars[0];
 			TestPositive(w0,"beam width");
-			vCopy(beam_pars+1,beam_center_0);
-			beam_asym=(beam_Npars==4 && (beam_center_0[0]!=0 || beam_center_0[1]!=0 || beam_center_0[2]!=0));
-			if (!beam_asym) vInit(beam_center);
 			s=1/(WaveNum*w0);
 			s2=s*s;
 			scale_x=1/w0;
@@ -181,10 +179,91 @@ void InitBeam(void)
 					default: LogError(ONE_POS,"Incompatibility error in GenerateB");
 				}
 				beam_descr=dyn_sprintf("Gaussian beam (%s)\n"
-				                       "\tWidth="GFORMDEF" (confinement factor s="GFORMDEF")\n"
-				                       "\tCenter position: "GFORMDEF3V,tmp_str,w0,s,COMP3V(beam_center_0));
+				                       "\tWidth="GFORMDEF" (confinement factor s="GFORMDEF")\n",tmp_str,w0,s);
 			}
 			return;
+#ifndef NO_FORTRAN
+		case B_BES_CS:
+		case B_BES_CSp:
+		case B_BES_M:
+		case B_BES_LE:
+		case B_BES_LM:
+		case B_BES_TEL:
+		case B_BES_TML:
+			if (surface) PrintError("Currently, Bessel incident beam is not supported for '-surf'");
+			// initialize parameters
+			ConvertToInteger(beam_pars[0],"beam order",&besN);
+			TestRangeII(besN,"beam order (might cause the incorrect calculation of Bessel function)",-50,50);
+			vorticity=besN;
+			besAlpha=Deg2Rad(beam_pars[1]);
+			besKt=WaveNum*sin(besAlpha);
+			besKz=WaveNum*cos(besAlpha);
+			switch (beamtype) { // definition of elements of matrix M ((Mex,Mey),(Mmx,Mmy))
+				case B_BES_CS:
+					TestRangeII(beam_pars[1],"half-cone angle",0,90);
+					besM[0]=0.5;  besM[1]=0;
+					besM[2]=0;    besM[3]=0.5;
+					if (IFROOT) tmp_str="circularly symmetric energy density";
+					break;
+				case B_BES_CSp:
+					TestRangeII(beam_pars[1],"half-cone angle",0,90);
+					besM[0]=0.5;  besM[1]=0;
+					besM[2]=0;    besM[3]=-0.5;
+					if (IFROOT) tmp_str="circularly symmetric energy density, alternative type";
+					break;
+				case B_BES_M:
+					TestRangeII(beam_pars[1],"half-cone angle",0,90);
+					if (beam_Npars==6) {
+						besM[0]=beam_pars[2];
+						besM[1]=beam_pars[3];
+						besM[2]=beam_pars[4];
+						besM[3]=beam_pars[5];
+					}
+					else {
+						besM[0]=beam_pars[2]+I*beam_pars[6];
+						besM[1]=beam_pars[3]+I*beam_pars[7];
+						besM[2]=beam_pars[4]+I*beam_pars[8];
+						besM[3]=beam_pars[5]+I*beam_pars[9];
+					}
+					if (IFROOT) tmp_str="generalized";
+					break;
+				case B_BES_LE:
+					TestRangeII(beam_pars[1],"half-cone angle",0,90);
+					besM[0]=0;  besM[1]=0;
+					besM[2]=0;  besM[3]=1;
+					if (IFROOT) tmp_str="linear electric field";
+					break;
+				case B_BES_LM:
+					TestRangeII(beam_pars[1],"half-cone angle",0,90);
+					besM[0]=0;  besM[1]=1;
+					besM[2]=0;  besM[3]=0;
+					if (IFROOT) tmp_str="linear magnetic field";
+					break;
+				/* TODO: for the following two types, both 0 and 90 degrees should be fine, but may require some
+				 * rearrangement of formulae. Then the tests for range of alpha should be moved outside of the case
+				 */
+				case B_BES_TEL:
+					TestRangeNN(beam_pars[1],"half-cone angle for TEL type",0,90);
+					besM[0]=-WaveNum/besKt;  besM[1]=0;
+					besM[2]=0;               besM[3]=besKz/besKt;
+					if (IFROOT) tmp_str="linear component of the TE";
+					break;
+				case B_BES_TML:
+					TestRangeNN(beam_pars[1],"half-cone angle for TML type",0,90);
+					besM[0]=0;              besM[1]=besKz/besKt;
+					besM[2]=WaveNum/besKt;  besM[2]=0;
+					if (IFROOT) tmp_str="linear component of the TM";
+					break;
+				default: LogError(ONE_POS,"Incompatibility error in GenerateB");
+			}
+			// TODO: some symmetries can be retained in some special cases
+			symR=symX=symY=symZ=false;
+			// beam info
+			if (IFROOT) beam_descr=dyn_sprintf("Bessel beam (%s)\n"
+				                               "\tOrder: %d, half-cone angle: "GFORMDEF" deg\n",
+				                               tmp_str,besN,beam_pars[1]);
+			return;
+#endif // !NO_FORTRAN
 		case B_READ:
 			// the safest is to assume cancellation of all symmetries
 			symX=symY=symZ=symR=false;
@@ -193,7 +272,7 @@ void InitBeam(void)
 				if (beam_Npars==1) beam_descr=dyn_sprintf("specified by file '%s'",beam_fnameY);
 				else beam_descr=dyn_sprintf("specified by files '%s' and '%s'",beam_fnameY,beam_fnameX);
 			}
-			// we do not define beam_asym here, because beam_center is not defined anyway
+			// this type is unaffected by beam_center
 			return;
 	}
 	LogError(ONE_POS,"Unknown type of incident beam (%d)",(int)beamtype);
@@ -208,11 +287,7 @@ void InitBeam(void)
 	 *    false here. Do not set any of them to true, as they can be set to false by other factors.
 	 *    symX, symY, symZ - symmetries of reflection over planes YZ, XZ, XY respectively.
 	 *    symR - symmetry of rotation for 90 degrees over the Z axis
-	 * 4) initialize the following:
-	 *    beam_descr - descriptive string, which will appear in log file.
-	 *    beam_asym - whether beam center does not coincide with the reference frame origin. If it is set to true, then
-	 *                set also beam_center_0 - 3D radius-vector of beam center in the laboratory reference frame (it
-	 *                will be then automatically transformed to particle reference frame, if required).
+	 * 4) initialize beam_descr - descriptive string, which will appear in log file.
 	 * 5) Consider the case of surface (substrate near the particle). If the new beam type is incompatible with it, add
 	 *    an explicit exception, like "if (surface) PrintErrorHelp(...);". Otherwise, you also need to define inc_scale.
 	 * All other auxiliary variables, which are used in beam generation (GenerateB(), see below), should be defined in
@@ -236,6 +311,21 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 	const double *ex; // coordinate axis of the beam reference frame
 	double ey[3];
 	double r1[3];
+	/* complex wave amplitudes of secondary waves (with phase relative to particle center);
+	 * The transmitted wave can be inhomogeneous wave (when msub is complex), then eIncTran (e) is normalized
+	 * counter-intuitively. Before multiplying by tc, it satisfies (e,e)=1!=||e||^2. This normalization is consistent
+	 * with used formulae for transmission coefficients. So this transmission coefficient is not (generally) equal to
+	 * the ratio of amplitudes of the electric fields. In particular, when E=E0*e, ||E||!=|E0|*||e||, where
+	 * ||e||^2=(e,e*)=|e_x|^2+|e_y|^2+|e_z|^2=1
+	 */
+	doublecomplex eIncRefl[3],eIncTran[3];
+#ifndef NO_FORTRAN
+	// for Bessel beams
+	int n1,q;
+	doublecomplex vort;  // vortex phase of Bessel beam rotated by 90 deg
+	doublecomplex fn[5]; // for general functions f(n,ro,phi) of Bessel beams (fn-2, fn-1, fn, fn+1, fn+2, respectively)
+	double phi,arg,td1[abs(besN)+3],td2[abs(besN)+3],jn1[abs(besN)+3]; // for Bessel beams
+#endif
 	const char *fname;
 	/* TO ADD NEW BEAM
 	 * Add here all intermediate variables, which are used only inside this function. You may as well use 't1'-'t8'
@@ -263,27 +353,23 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 				 * original medium.
 				 */
 				doublecomplex rc,tc; // reflection and transmission coefficients
+				double hbeam=hsub+beam_center[2]; // height of beam center above the surface 
 				if (prop[2]>0) { // beam comes from the substrate (below)
 					//  determine amplitude of the reflected and transmitted waves; here msub is always defined
 					if (which==INCPOL_Y) { // s-polarized
-						cvBuildRe(ex,eIncRefl);
 						cvBuildRe(ex,eIncTran);
-						rc=FresnelRS(ki,kt);
 						tc=FresnelTS(ki,kt);
 					}
 					else { // p-polarized
-						vInvRefl_cr(ex,eIncRefl);
 						crCrossProd(ey,ktVec,eIncTran);
-						rc=FresnelRP(ki,kt,1/msub);
 						tc=FresnelTP(ki,kt,1/msub);
 					}
-					// phase shift due to the origin at height hsub
-					cvMultScal_cmplx(rc*cexp(-2*I*WaveNum*ki*hsub),eIncRefl,eIncRefl);
-					cvMultScal_cmplx(tc*cexp(I*WaveNum*(kt-ki)*hsub),eIncTran,eIncTran);
+					// phase shift due to the beam center relative to the origin and surface
+					cvMultScal_cmplx(tc*cexp(I*WaveNum*((kt-ki)*hbeam-crDotProd(ktVec,beam_center))),eIncTran,eIncTran);
 					// main part
 					for (i=0;i<local_nvoid_Ndip;i++) {
 						j=3*i;
-						// b[i] = eIncTran*exp(ik*kt.r)
+						// b[i] = eIncTran*exp(ik*kt.r); subtraction of beam_center is avoided by the factor above
 						cvMultScal_cmplx(cexp(I*WaveNum*crDotProd(ktVec,DipoleCoord+j)),eIncTran,b+j);
 					}
 				}
@@ -291,44 +377,38 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 					// determine amplitude of the reflected and transmitted waves
 					if (which==INCPOL_Y) { // s-polarized
 						cvBuildRe(ex,eIncRefl);
-						if (msubInf) {
-							rc=-1;
-							tc=0; // to remove compiler warnings
-						}
+						if (msubInf) rc=-1;
 						else {
 							cvBuildRe(ex,eIncTran);
 							rc=FresnelRS(ki,kt);
-							tc=FresnelTS(ki,kt);
 						}
 					}
 					else { // p-polarized
 						vInvRefl_cr(ex,eIncRefl);
-						if (msubInf) {
-							rc=1;
-							tc=0; // to remove compiler warnings
-						}
+						if (msubInf) rc=1;
 						else {
 							crCrossProd(ey,ktVec,eIncTran);
 							cvMultScal_cmplx(1/msub,eIncTran,eIncTran); // normalize eIncTran by ||ktVec||=msub
 							rc=FresnelRP(ki,kt,msub);
-							tc=FresnelTP(ki,kt,msub);
 						}
 					}
-					// phase shift due to the origin at height hsub
-					cvMultScal_cmplx(rc*imExp(2*WaveNum*creal(ki)*hsub),eIncRefl,eIncRefl); // assumes real ki
-					if (!msubInf) cvMultScal_cmplx(tc*cexp(I*WaveNum*(ki-kt)*hsub),eIncTran,eIncTran);
+					// phase shift due to the beam center relative to the surface
+					cvMultScal_cmplx(rc*imExp(2*WaveNum*creal(ki)*hbeam),eIncRefl,eIncRefl); // assumes real ki
 					// main part
 					for (i=0;i<local_nvoid_Ndip;i++) {
 						j=3*i;
+						// together with above phase shift, beam_center affects only the common phase factor
+						LinComb(DipoleCoord+j,beam_center,1,-1,r1); 
 						// b[i] = ex*exp(ik*r.a) + eIncRefl*exp(ik*prIncRefl.r)
-						cvMultScal_RVec(imExp(WaveNum*DotProd(DipoleCoord+j,prop)),ex,b+j);
-						cvLinComb1_cmplx(eIncRefl,b+j,imExp(WaveNum*DotProd(DipoleCoord+j,prIncRefl)),b+j);
+						cvMultScal_RVec(imExp(WaveNum*DotProd(r1,prop)),ex,b+j);
+						cvLinComb1_cmplx(eIncRefl,b+j,imExp(WaveNum*DotProd(r1,prIncRefl)),b+j);
 					}
 				}
 			}
 			else for (i=0;i<local_nvoid_Ndip;i++) { // standard (non-surface) plane wave
 				j=3*i;
-				ctemp=imExp(WaveNum*DotProd(DipoleCoord+j,prop)); // ctemp=exp(ik*r.a)
+				LinComb(DipoleCoord+j,beam_center,1,-1,r1); // this affects only the common phase factor
+				ctemp=imExp(WaveNum*DotProd(r1,prop)); // ctemp=exp(ik*r.a)
 				cvMultScal_RVec(ctemp,ex,b+j); // b[i]=ctemp*ex
 			}
 			return;
@@ -438,6 +518,105 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 				}
 			}
 			return;
+#ifndef NO_FORTRAN
+		case B_BES_CS:
+		case B_BES_CSp:
+		case B_BES_M:
+		case B_BES_LE:
+		case B_BES_LM:
+		case B_BES_TEL:
+		case B_BES_TML:
+			/* we assume that matrix M determines x-polarization (in the beam reference frame), while y-one is obtained
+			 * by rotation with additional phase shift
+			 */
+			vort=(which==INCPOL_Y) ? cpow(I,besN) : 1;
+			for (i=0;i<local_nvoid_Ndip;i++) {
+				j=3*i;
+				LinComb(DipoleCoord+j,beam_center,1,-1,r1);
+				x=DotProd(r1,ex);
+				y=DotProd(r1,ey);
+				z=DotProd(r1,prop);
+				phi=atan2(y,x); // angular coordinate in a cylindrical coordinate system
+				ctemp=imExp(besN*phi)*cexp(I*besKz*z)*vort/(WaveNum*WaveNum); // common factor
+				arg=besKt*sqrt(x*x+y*y); // argument of Bessel functions
+				if (arg<ROUND_ERR) {
+					// TODO: the following seems incorrect for n=+-1 or +-2 (other fn will be non-zero)
+					if (besN==0) fn[2]=1.;
+					else fn[2]=0;
+					fn[0]=fn[1]=fn[3]=fn[4]=0;
+				}
+				else {
+					// TODO: the following looks very complicated, try to simplify
+					n1=abs(besN)+2;
+					if (besN<=-3) {
+						bjndd_(&n1,&arg,jn1,td1,td2);
+						if (n1%2==0) q=1;
+						else q=-1;
+						fn[0] =  q*jn1[-besN+2]*imExp(-2*phi);
+						fn[1] = -q*jn1[-besN+1]*imExp(-phi);
+						fn[2] =  q*jn1[-besN];
+						fn[3] = -q*jn1[-besN-1]*imExp(phi);
+						fn[4] =  q*jn1[-besN-2]*imExp(2*phi);
+					}
+					if (besN >= 2) {
+						bjndd_(&n1,&arg,jn1,td1,td2);
+						fn[0] = jn1[besN-2]*imExp(-2*phi);
+						fn[1] = jn1[besN-1]*imExp(-phi);
+						fn[2] = jn1[besN];
+						fn[3] = jn1[besN+1]*imExp(phi);
+						fn[4] = jn1[besN+2]*imExp(2*phi);
+					}
+					if (besN == -2) {
+						bjndd_(&n1,&arg,jn1,td1,td2);
+						fn[0] =  jn1[4]*imExp(-2*phi);
+						fn[1] = -jn1[3]*imExp(-phi);
+						fn[2] =  jn1[2];
+						fn[3] = -jn1[1]*imExp(phi);
+						fn[4] =  jn1[0]*imExp(2*phi);
+					}
+					if (besN == -1) {
+						bjndd_(&n1,&arg,jn1,td1,td2);
+						fn[0] = -jn1[3]*imExp(-2*phi);
+						fn[1] =  jn1[2]*imExp(-phi);
+						fn[2] = -jn1[1];
+						fn[3] =  jn1[0]*imExp(phi);
+						fn[4] =  jn1[1]*imExp(2*phi);
+					}
+					if (besN == 0) {
+						bjndd_(&n1,&arg,jn1,td1,td2);
+						fn[0] =  jn1[2]*imExp(-2*phi);
+						fn[1] = -jn1[1]*imExp(-phi);
+						fn[2] =  jn1[0];
+						fn[3] =  jn1[1]*imExp(phi);
+						fn[4] =  jn1[2]*imExp(2*phi);
+					}
+					if (besN == 1) {
+						bjndd_(&n1,&arg,jn1,td1,td2);
+						fn[0] = -jn1[1]*imExp(-2*phi);
+						fn[1] =  jn1[0]*imExp(-phi);
+						fn[2] =  jn1[1];
+						fn[3] =  jn1[2]*imExp(phi);
+						fn[4] =  jn1[3]*imExp(2*phi);
+					}
+				}
+				/* TODO: all factors before f[n] should be calculated in InitBeam beforehand
+				 * this will improve speed and allow robust calculation in the limit, e.g., of kt->0
+				 */
+				t1 = (((WaveNum*WaveNum+besKz*besKz)/2.*besM[0] +  WaveNum*besKz*besM[3])*fn[2] +
+					  besKt*besKt/4.*((besM[0]+I*besM[1])*fn[0] + (besM[0]-I*besM[1])*fn[4])); // Ex
+				t2 = (((WaveNum*WaveNum+besKz*besKz)/2.*besM[1] -  WaveNum*besKz*besM[2])*fn[2] +
+					  I*besKt*besKt/4.*((besM[0]+I*besM[1])*fn[0] - (besM[0]-I*besM[1])*fn[4])); // Ey
+				t3 = ((I*besKz*(besM[0]+I*besM[1]) + WaveNum*(besM[2]+I*besM[3]))*fn[1] -
+					  (I*besKz*(besM[0]-I*besM[1]) - WaveNum*(besM[2]-I*besM[3]))*fn[3])*besKt/2.; // Ez
+
+				cvMultScal_RVec(t1,ex,v1);
+				cvMultScal_RVec(t2,ey,v2);
+				cvMultScal_RVec(t3,prop,v3);
+				cvAdd2Self(v1,v2,v3);
+				cvMultScal_cmplx(ctemp,v1,b+j);
+			}
+			return;
+#endif // !NO_FORTRAN		
 		case B_READ:
 			if (which==INCPOL_Y) fname=beam_fnameY;
 			else fname=beam_fnameX; // which==INCPOL_X
@@ -449,12 +628,12 @@ void GenerateB (const enum incpol which,   // x - or y polarized incident light
 	 * add a case above. Identifier ('B_...') should be defined inside 'enum beam' in const.h. This case should set
 	 * complex vector 'b', describing the incident field in the particle reference frame. It is set inside the cycle for
 	 * each dipole of the particle and is calculated using
-	 * 1) 'DipoleCoord' � array of dipole coordinates;
-	 * 2) 'prop' � propagation direction of the incident field;
-	 * 3) 'ex' � direction of incident polarization;
-	 * 4) 'ey' � complementary unity vector of polarization (orthogonal to both 'prop' and 'ex');
-	 * 5) 'beam_center' � beam center in the particle reference frame (automatically calculated from 'beam_center_0'
-	 *                    defined in InitBeam).
+	 * 1) 'DipoleCoord' - array of dipole coordinates;
+	 * 2) 'prop' - propagation direction of the incident field;
+	 * 3) 'ex' - direction of incident polarization;
+	 * 4) 'ey' - complementary unity vector of polarization (orthogonal to both 'prop' and 'ex');
+	 * 5) 'beam_center' - beam center in the particle reference frame (automatically calculated from 'beam_center_0'
+	 *                    defined by '-beam_center' command line option).
 	 * If the new beam type is compatible with '-surf', include here the corresponding code. For that you will need
 	 * the variables, related to surface - see vars.c after "// related to a nearby surface".
 	 * If you need temporary local variables (which are used only in this part of the code), either use 't1'-'t8' or
diff --git a/src/Makefile b/src/Makefile
index 0380a433..f23b53bd 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -192,7 +192,7 @@ CSOURCE := ADDAmain.c CalculateE.c calculator.c chebyshev.c cmplx.c comm.c cross
 # Fortran files are located in src/fort folder, other files may be added below
 FSOURCE := d07hre.f d09hre.f d113re.f d132re.f dadhre.f dchhre.f dcuhre.f dfshre.f dinhre.f drlhre.f dtrhre.f \
            propaesplibreintadda.f
-F90SOURCE := 
+F90SOURCE := bessel.f90
 # C++ files are located in src/cpp folder, other files may be added below
 CPPSOURCE :=
 
diff --git a/src/README.md b/src/README.md
index 74437518..4d32ea25 100644
--- a/src/README.md
+++ b/src/README.md
@@ -1,6 +1,7 @@
 Main source of ADDA, including Makefiles and other scripts. C++ and Fortran sources of auxiliary (optional) routines are placed in separate folders. The compilation itself happens in `seq`, `mpi`, or `ocl` folder depending on the mode - object and executable files are placed there.
 * `cpp/` - C++ sources for Apple clFFT
 * `fort/` - Fortran sources
+  * `bessel.f90` - subroutines for calculating Bessel functions
   * `cfft99D.f` - source file for Temperton FFT
   * `propaesplibreintadda.f`, `d07hre.f`, `d09hre.f`, `d113re.f`, `d132re.f`, `dadhre.f`, `dchhre.f`, `dcuhre.f`, `dfshre.f`, `dinhre.f`, `drlhre.f`, `dtrhre.f` - subroutine to compute IGT and related numerical routines
 * `mpi/Makefile` - makefile for MPI version (called from the main makefile)
diff --git a/src/const.h b/src/const.h
index 5d655710..b0748e48 100644
--- a/src/const.h
+++ b/src/const.h
@@ -285,6 +285,15 @@ enum incpol {
 
 enum beam { // beam types
 	B_BARTON5, // 5th order description of the Gaussian beam
+#ifndef NO_FORTRAN
+	B_BES_CS,  // Bessel beam with circularly symmetric energy density
+	B_BES_CSp, // Alternative Bessel beam with circularly symmetric energy density
+	B_BES_M,   // Generalized Bessel beam
+	B_BES_LE,  // Bessel beam with linearly polarized electric field
+	B_BES_LM,  // Bessel beam with linearly polarized magnetic field
+	B_BES_TEL, // Linear component of the TE Bessel beam
+	B_BES_TML, // Linear component of the TM Bessel beam
+#endif
 	B_DAVIS3,  // 3rd order description of the Gaussian beam
 	B_DIPOLE,  // field of a point dipole
 	B_LMINUS,  // 1st order description of the Gaussian beam
diff --git a/src/crosssec.c b/src/crosssec.c
index c62b23fb..31499acc 100644
--- a/src/crosssec.c
+++ b/src/crosssec.c
@@ -40,9 +40,6 @@ extern const doublecomplex cc[][3];
 #ifndef SPARSE
 extern doublecomplex * restrict expsX,* restrict expsY,* restrict expsZ;
 #endif
-// defined and initialized in GenerateB.c
-extern const double beam_center_0[3];
-//extern doublecomplex eIncRefl[3],eIncTran[3];
 // defined and initialized in param.c
 extern const double incPolX_0[3],incPolY_0[3];
 extern const enum scat ScatRelation;
@@ -831,7 +828,8 @@ double ExtCross(const double * restrict incPol)
 	double sum;
 	size_t i;
 
-	if (beamtype==B_PLANE && !surface) {
+	// this can be considered a legacy case, which works only for the simplest plane way centered at the particle 
+	if (beamtype==B_PLANE && !surface && !beam_asym) {
 		CalcField (ebuff,prop);
 		sum=crDotProd_Re(ebuff,incPol); // incPol is real, so no conjugate is needed
 		MyInnerProduct(&sum,double_type,1,&Timing_ScatQuanComm);
@@ -850,8 +848,8 @@ double ExtCross(const double * restrict incPol)
 		 * independent of propagation or scattering direction. For rectangular dipoles, it is only approximate but
 		 * expected to be accurate for not very elongated dipoles and/or not very spread out incident field.
 		 *
-		 * In principle, the situation is similar for SO of full IGT, but there the correction factor depends on the
-		 * propagation direction.even for cubical dipoles
+		 * In principle, the situation is similar for full IGT, but there the correction factor depends on the
+		 * propagation direction even for cubical dipoles
 		 */
 		if (ScatRelation==SQ_IGT_SO) sum*=eta2(prop);
 	}
diff --git a/src/fort/bessel.f90 b/src/fort/bessel.f90
new file mode 100644
index 00000000..6609b381
--- /dev/null
+++ b/src/fort/bessel.f90
@@ -0,0 +1,304 @@
+subroutine bjndd ( n, x, bj, dj, fj )
+
+!*****************************************************************************80
+!
+!! BJNDD computes Bessel functions Jn(x) and first and second derivatives.
+!
+!  Licensing:
+!
+!    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
+!    they give permission to incorporate this routine into a user program
+!    provided that the copyright is acknowledged.
+!
+!  Modified:
+!
+!    11 July 2012
+!
+!  Author:
+!
+!    Shanjie Zhang, Jianming Jin
+!
+!  Reference:
+!
+!    Shanjie Zhang, Jianming Jin,
+!    Computation of Special Functions,
+!    Wiley, 1996,
+!    ISBN: 0-471-11963-6,
+!    LC: QA351.C45.
+!
+!  Parameters:
+!
+!    Input, integer ( kind = 4 ) N, the order.
+!
+!    Input, real ( kind = 8 ) X, the argument.
+!
+!    Output, real ( kind = 8 ) BJ(N+1), DJ(N+1), FJ(N+1), the values of
+!    Jn(x), Jn'(x) and Jn''(x) in the last entries.
+!
+  implicit none
+
+  integer ( kind = 4 ) n
+
+  real ( kind = 8 ) bj(n+1)
+  real ( kind = 8 ) bs
+  real ( kind = 8 ) dj(n+1)
+  real ( kind = 8 ) f
+  real ( kind = 8 ) f0
+  real ( kind = 8 ) f1
+  real ( kind = 8 ) fj(n+1)
+  integer ( kind = 4 ) k
+  integer ( kind = 4 ) m
+  integer ( kind = 4 ) mt
+  integer ( kind = 4 ) nt
+  real ( kind = 8 ) x
+
+  do nt = 1, 900
+    mt = int ( 0.5D+00 * log10 ( 6.28D+00 * nt ) &
+      - nt * log10 ( 1.36D+00 * abs ( x ) / nt ) )
+    if ( 20 < mt ) then
+      exit
+    end if
+  end do
+
+  m = nt
+  bs = 0.0D+00
+  f0 = 0.0D+00
+  f1 = 1.0D-35
+  do k = m, 0, -1
+    f = 2.0D+00 * ( k + 1.0D+00 ) * f1 / x - f0
+    if ( k <= n ) then
+      bj(k+1) = f
+    end if
+    if ( k == 2 * int ( k / 2 ) ) then
+      bs = bs + 2.0D+00 * f
+    end if
+    f0 = f1
+    f1 = f
+  end do
+
+  do k = 0, n
+    bj(k+1) = bj(k+1) / ( bs - f )
+  end do
+
+  dj(1) = -bj(2)
+  fj(1) = -1.0D+00 * bj(1) - dj(1) / x
+  do k = 1, n
+    dj(k+1) = bj(k) - k * bj(k+1) / x
+    fj(k+1) = ( k * k / ( x * x ) - 1.0D+00 ) * bj(k+1) - dj(k+1) / x
+  end do
+
+  return
+end
+subroutine cik01 ( z, cbi0, cdi0, cbi1, cdi1, cbk0, cdk0, cbk1, cdk1 )
+
+!*****************************************************************************80
+!
+!! CIK01: modified Bessel I0(z), I1(z), K0(z) and K1(z) for complex argument.
+!
+!  Discussion:
+!
+!    This procedure computes the modified Bessel functions I0(z), I1(z),
+!    K0(z), K1(z), and their derivatives for a complex argument.
+!
+!  Licensing:
+!
+!    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
+!    they give permission to incorporate this routine into a user program
+!    provided that the copyright is acknowledged.
+!
+!  Modified:
+!
+!    31 July 2012
+!
+!  Author:
+!
+!    Shanjie Zhang, Jianming Jin
+!
+!  Reference:
+!
+!    Shanjie Zhang, Jianming Jin,
+!    Computation of Special Functions,
+!    Wiley, 1996,
+!    ISBN: 0-471-11963-6,
+!    LC: QA351.C45.
+!
+!  Parameters:
+!
+!    Input, complex ( kind = 8 ) Z, the argument.
+!
+!    Output, complex ( kind = 8 ) CBI0, CDI0, CBI1, CDI1, CBK0, CDK0, CBK1,
+!    CDK1, the values of I0(z), I0'(z), I1(z), I1'(z), K0(z), K0'(z), K1(z),
+!    and K1'(z).
+!
+  implicit none
+
+  real ( kind = 8 ), save, dimension ( 12 ) :: a = (/ &
+    0.125D+00,           7.03125D-02,&
+    7.32421875D-02,      1.1215209960938D-01,&
+    2.2710800170898D-01, 5.7250142097473D-01,&
+    1.7277275025845D+00, 6.0740420012735D+00,&
+    2.4380529699556D+01, 1.1001714026925D+02,&
+    5.5133589612202D+02, 3.0380905109224D+03 /)
+  real ( kind = 8 ) a0
+  real ( kind = 8 ), save, dimension ( 10 ) :: a1 = (/ &
+    0.125D+00,            0.2109375D+00, &
+    1.0986328125D+00,     1.1775970458984D+01, &
+    2.1461706161499D+002, 5.9511522710323D+03, &
+    2.3347645606175D+05,  1.2312234987631D+07, &
+    8.401390346421D+08,   7.2031420482627D+10 /)
+  real ( kind = 8 ), save, dimension ( 12 ) :: b = (/ &
+   -0.375D+00,           -1.171875D-01, &
+   -1.025390625D-01,     -1.4419555664063D-01, &
+   -2.7757644653320D-01, -6.7659258842468D-01, &
+   -1.9935317337513D+00, -6.8839142681099D+00, &
+   -2.7248827311269D+01, -1.2159789187654D+02, &
+   -6.0384407670507D+02, -3.3022722944809D+03 /)
+  complex ( kind = 8 ) ca
+  complex ( kind = 8 ) cb
+  complex ( kind = 8 ) cbi0
+  complex ( kind = 8 ) cbi1
+  complex ( kind = 8 ) cbk0
+  complex ( kind = 8 ) cbk1
+  complex ( kind = 8 ) cdi0
+  complex ( kind = 8 ) cdi1
+  complex ( kind = 8 ) cdk0
+  complex ( kind = 8 ) cdk1
+  complex ( kind = 8 ) ci
+  complex ( kind = 8 ) cr
+  complex ( kind = 8 ) cs
+  complex ( kind = 8 ) ct
+  complex ( kind = 8 ) cw
+  integer ( kind = 4 ) k
+  integer ( kind = 4 ) k0
+  real ( kind = 8 ) pi
+  real ( kind = 8 ) w0
+  complex ( kind = 8 ) z
+  complex ( kind = 8 ) z1
+  complex ( kind = 8 ) z2
+  complex ( kind = 8 ) zr
+  complex ( kind = 8 ) zr2
+
+  pi = 3.141592653589793D+00
+  ci = cmplx ( 0.0D+00, 1.0D+00, kind = 8 )
+  a0 = abs ( z )
+  z2 = z * z
+  z1 = z
+
+  if ( a0 == 0.0D+00 ) then
+    cbi0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    cbi1 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
+    cdi0 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
+    cdi1 = cmplx ( 0.5D+00, 0.0D+00, kind = 8 )
+    cbk0 = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
+    cbk1 = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
+    cdk0 = - cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
+    cdk1 = - cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
+    return
+  end if
+
+  if ( real ( z, kind = 8 ) < 0.0D+00 ) then
+    z1 = -z
+  end if
+
+  if ( a0 <= 18.0D+00 ) then
+
+    cbi0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    do k = 1, 50
+      cr = 0.25D+00 * cr * z2 / ( k * k )
+      cbi0 = cbi0 + cr
+      if ( abs ( cr / cbi0 ) < 1.0D-15 ) then
+        exit
+      end if
+    end do
+
+    cbi1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    do k = 1, 50
+      cr = 0.25D+00 * cr * z2 / ( k * ( k + 1 ) )
+      cbi1 = cbi1 + cr
+      if ( abs ( cr / cbi1 ) < 1.0D-15 ) then
+        exit
+      end if
+    end do
+
+    cbi1 = 0.5D+00 * z1 * cbi1
+
+  else
+
+    if ( a0 < 35.0D+00 ) then
+      k0 = 12
+    else if ( a0 < 50.0D+00 ) then
+      k0 = 9
+    else
+      k0 = 7
+    end if
+
+    ca = exp ( z1 ) / sqrt ( 2.0D+00 * pi * z1 )
+    cbi0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    zr = 1.0D+00 / z1
+    do k = 1, k0
+      cbi0 = cbi0 + a(k) * zr ** k
+    end do
+    cbi0 = ca * cbi0
+    cbi1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    do k = 1, k0
+      cbi1 = cbi1 + b(k) * zr ** k
+    end do
+    cbi1 = ca * cbi1
+
+  end if
+
+  if ( a0 <= 9.0D+00 ) then
+
+    cs = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
+    ct = - log ( 0.5D+00 * z1 ) - 0.5772156649015329D+00
+    w0 = 0.0D+00
+    cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    do k = 1, 50
+      w0 = w0 + 1.0D+00 / k
+      cr = 0.25D+00 * cr / ( k * k ) * z2
+      cs = cs + cr * ( w0 + ct )
+      if ( abs ( ( cs - cw ) / cs ) < 1.0D-15 ) then
+        exit
+      end if
+      cw = cs
+    end do
+
+    cbk0 = ct + cs
+
+  else
+
+    cb = 0.5D+00 / z1
+    zr2 = 1.0D+00 / z2
+    cbk0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
+    do k = 1, 10
+      cbk0 = cbk0 + a1(k) * zr2 ** k
+    end do
+    cbk0 = cb * cbk0 / cbi0
+
+  end if
+
+  cbk1 = ( 1.0D+00 / z1 - cbi1 * cbk0 ) / cbi0
+
+  if ( real ( z, kind = 8 ) < 0.0D+00 ) then
+
+    if ( imag ( z ) < 0.0D+00 ) then
+      cbk0 = cbk0 + ci * pi * cbi0
+      cbk1 = - cbk1 + ci * pi * cbi1
+    else
+      cbk0 = cbk0 - ci * pi * cbi0
+      cbk1 = - cbk1 - ci * pi * cbi1
+    end if
+
+    cbi1 = - cbi1
+
+  end if
+
+  cdi0 = cbi1
+  cdi1 = cbi0 - 1.0D+00 / z * cbi1
+  cdk0 = - cbk1
+  cdk1 = - cbk0 - 1.0D+00 / z * cbk1
+
+  return
+end
diff --git a/src/param.c b/src/param.c
index 12582533..e5767bac 100644
--- a/src/param.c
+++ b/src/param.c
@@ -188,6 +188,8 @@ static bool orient_used;        // whether '-orient ...' was used in the command
 static bool yz_used;            // whether '-yz ...' was used in the command line
 static bool scat_plane_used;    // whether '-scat_plane ...' was used in the command line
 static bool so_buf_used;        // whether '-so_buf ...' was used in the command line
+static bool beam_center_used;   // whether '-beam_center ...' was used in the command line
+static bool deprecated_bc_used; // whether '-beam ... <x> <y> <z>' was used in the command line (deprecated option)
 
 /* TO ADD NEW COMMAND LINE OPTION
  * If you need new variables or flags to implement effect of the new command line option, define them here. If a
@@ -228,17 +230,38 @@ static const char exeusage[]="[-<opt1> [<args1>] [-<opt2> <args2>]...]]";
 static const struct subopt_struct beam_opt[]={
 	{"barton5","<width> [<x> <y> <z>]","5th order approximation of the Gaussian beam (by Barton). The beam width is "
 		"obligatory and x, y, z coordinates of the center of the beam (in laboratory reference frame) are optional "
-		"(zero, by default). All arguments are in um. This is recommended option for simulation of the Gaussian beam.",
-		UNDEF,B_BARTON5},
+		"(zero, by default). All arguments are in um. This is recommended option for simulation of the Gaussian beam. "
+		"Specification of coordinates here is DEPRECATED, use -beam_center instead.",UNDEF,B_BARTON5},
+#ifndef NO_FORTRAN
+	{"besselCS","<order> <angle>","Bessel beam with circularly symmetric energy density. Order is integer (of any "
+		"sign) and the half-cone angle (in degrees) is measured from the z-axis.",2,B_BES_CS},
+	{"besselCSp","<order> <angle>","Alternative Bessel beam with circularly symmetric energy density. Order is "
+		"integer (of any sign) and the half-cone angle (in degrees) is measured from the z-axis.",2,B_BES_CSp},
+	{"besselM","<order> <angle> <ReMex> <ReMey> <ReMmx> <ReMmy> [<ImMex> <ImMey> <ImMmx> <ImMmy>]",
+		"Generalized Bessel beam. Order is integer (of any sign) and the half-cone angle (in degrees) is measured from "
+		"the z-axis. The beam is defined by 2x2 matrix M: (Mex, Mey, Mmx, Mmy). Real parts of these four elements are "
+		"obligatory, while imaginary parts are optional (zero, by default).",UNDEF,B_BES_M},
+	{"besselLE","<order> <angle>","Bessel beam with linearly polarized electric field. Order is integer (of any sign) "
+		"and the half-cone angle (in degrees) is measured from the z-axis.",2,B_BES_LE},
+	{"besselLM","<order> <angle>","Bessel beam with linearly polarized magnetic field. Order is integer (of any sign) "
+		"and the half-cone angle (in degrees) is measured from the z-axis.",2,B_BES_LM},
+	{"besselTEL","<order> <angle>","Linear component of the TE Bessel beam. Order is integer (of any sign) and the "
+		"half-cone angle (in degrees) is measured from the z-axis.",2,B_BES_TEL},
+	{"besselTML","<order> <angle>","Linear component of the TM Bessel beam. Order is integer (of any sign) and the "
+		"half-cone angle (in degrees) is measured from the z-axis.",2,B_BES_TML},
+#endif // !NO_FORTRAN
 	{"davis3","<width> [<x> <y> <z>]","3rd order approximation of the Gaussian beam (by Davis). The beam width is "
 		"obligatory and x, y, z coordinates of the center of the beam (in laboratory reference frame) are optional "
-		"(zero, by default). All arguments are in um.",UNDEF,B_DAVIS3},
-	{"dipole","<x> <y> <z>","Field of a unit point dipole placed at x, y, z coordinates (in laboratory reference "
+		"(zero, by default). All arguments are in um. Specification of coordinates here is DEPRECATED, use "
+		"-beam_center instead.",UNDEF,B_DAVIS3},
+	{"dipole","[<x> <y> <z>]","Field of a unit point dipole placed at x, y, z coordinates (in laboratory reference "
 		"frame). All arguments are in um. Orientation of the dipole is determined by -prop command line option."
-		"Implies '-scat_matr none'. If '-surf' is used, dipole position should be above the surface.",3,B_DIPOLE},
+		"Implies '-scat_matr none'. If '-surf' is used, dipole position should be above the surface. Specification of "
+		"coordinates here is DEPRECATED, use -beam_center instead.",UNDEF,B_DIPOLE},
 	{"lminus","<width> [<x> <y> <z>]","Simplest approximation of the Gaussian beam. The beam width is obligatory and "
 		"x, y, z coordinates of the center of the beam (in laboratory reference frame) are optional (zero, by"
-		" default). All arguments are in um.",UNDEF,B_LMINUS},
+		" default). All arguments are in um. Specification of coordinates here is DEPRECATED, use -beam_center "
+		"instead.",UNDEF,B_LMINUS},
 	{"plane","","Infinite plane wave",0,B_PLANE},
 	{"read","<filenameY> [<filenameX>]","Defined by separate files, which names are given as arguments. Normally two "
 		"files are required for Y- and X-polarizations respectively, but a single filename is sufficient if only "
@@ -339,6 +362,7 @@ PARSE_FUNC(alldir_inp);
 PARSE_FUNC(anisotr);
 PARSE_FUNC(asym);
 PARSE_FUNC(beam);
+PARSE_FUNC(beam_center);
 PARSE_FUNC(chp_dir);
 PARSE_FUNC(chp_load);
 PARSE_FUNC(chp_type);
@@ -417,8 +441,12 @@ static struct opt_struct options[]={
 		"formulations, and '-rect_dip'.",0,NULL},
 	{PAR(asym),"","Calculate the asymmetry vector. Implies '-Csca' and '-vec'",0,NULL},
 	{PAR(beam),"<type> [<args>]","Sets the incident beam, either predefined or 'read' from file. All parameters of "
-		"predefined beam types (if present) are floats.\n"
+		"predefined beam types are floats except for <order> or filenames.\n"
 		"Default: plane",UNDEF,beam_opt},
+	{PAR(beam_center),"<x> <y> <z>","Sets the center of the beam in the laboratory reference frame (in um). For most "
+		"beams it corresponds to the most symmetric point with zero phase, while for a point source or a fast "
+		"electron, it determines the real position in space.\n"
+		"Default: 0 0 0",3,NULL},
 	{PAR(chp_dir),"<dirname>","Sets directory for the checkpoint (both for saving and loading).\n"
 		"Default: "FD_CHP_DIR,1,NULL},
 	{PAR(chp_load),"","Restart a simulation from a checkpoint",0,NULL},
@@ -634,7 +662,7 @@ static struct opt_struct options[]={
 		 * the next string.
 		 */
 	{PAR(shape),"<type> [<args>]","Sets shape of the particle, either predefined or 'read' from file. All parameters "
-		"of predefined shapes are floats except for filenames.\n"
+		"of predefined shapes are floats except for <n> and filenames.\n"
 		"Default: sphere",UNDEF,shape_opt},
 	{PAR(size),"<arg>","Sets the size of the computational grid along the x-axis in um, float. If default wavelength "
 		"is used, this option specifies the 'size parameter' of the computational grid. Can not be used together with "
@@ -856,6 +884,16 @@ static void ScanDoubleError(const char * restrict str,double *res)
 
 //======================================================================================================================
 
+static void ScanDouble3Error(char **argv,double res[static 3])
+// scans an option argument (3D vector of doubles) and checks for errors
+{
+	ScanDoubleError(argv[0],res);
+	ScanDoubleError(argv[1],res+1);
+	ScanDoubleError(argv[2],res+2);
+}
+
+//======================================================================================================================
+
 static void ScanIntError(const char * restrict str,int *res)
 // scanf an option argument and checks for errors
 {
@@ -994,11 +1032,21 @@ PARSE_FUNC(beam)
 		beam_Npars=Narg;
 		opt_beam=opt;
 		need=beam_opt[i].narg;
-		// check number of arguments
+		// check number of arguments and process deprecated ones
 		switch (beamtype) {
 			case B_LMINUS:
 			case B_DAVIS3:
-			case B_BARTON5: if (Narg!=1 && Narg!=4) NargError(Narg,"1 or 4"); break;
+			case B_BARTON5:
+				if (Narg!=1 && Narg!=4) NargError(Narg,"1 or 4");
+				if (Narg==4) deprecated_bc_used=true;
+				break;
+			case B_DIPOLE:
+				if (Narg!=0 && Narg!=3) NargError(Narg,"0 or 3");
+				if (Narg==3) deprecated_bc_used=true;
+				break;
+#ifndef NO_FORTRAN
+			case B_BES_M: if (Narg!=6 && Narg!=10) NargError(Narg,"6 or 10"); break;
+#endif
 			default: TestNarg(Narg,need); break;
 		}
 		/* TO ADD NEW BEAM
@@ -1010,10 +1058,20 @@ PARSE_FUNC(beam)
 			for (j=0;j<Narg;j++) ScanDoubleError(argv[j+2],beam_pars+j);
 		// stop search
 		found=true;
+		if (deprecated_bc_used) {
+			LogWarning(EC_WARN,ONE_POS,
+				"Providing beam-center coordinates as arguments to '-beam' is deprecated. Use '-beam_center' instead.");
+			vCopy(beam_pars+Narg-3,beam_center_0);
+		}
 		break;
 	}
 	if(!found) NotSupported("Beam type",argv[1]);
 }
+PARSE_FUNC(beam_center)
+{
+	ScanDouble3Error(argv+1,beam_center_0);
+	beam_center_used = true;
+}
 PARSE_FUNC(chp_dir)
 {
 	chp_dir=ScanStrError(argv[1],MAX_DIRNAME);
@@ -1140,6 +1198,10 @@ PARSE_FUNC(h)
 #ifdef SPARSE
 						if (strcmp(options[i].name,"shape")==0)
 							printf("!!! Most of the shape options are disabled in sparse mode\n");
+#endif
+#ifdef NO_FORTRAN
+						if (strcmp(options[i].name,"beam")==0)
+							printf("!!! Bessel beams rely on Fortran sources that were disabled at compile time.\n");
 #endif
 					}
 				}
@@ -1395,9 +1457,7 @@ PARSE_FUNC(prop)
 {
 	double tmp;
 
-	ScanDoubleError(argv[1],prop_0);
-	ScanDoubleError(argv[2],prop_0+1);
-	ScanDoubleError(argv[3],prop_0+2);
+	ScanDouble3Error(argv+1,prop_0);
 	tmp=DotProd(prop_0,prop_0);
 	if (tmp==0) PrintErrorHelp("Given propagation vector is null");
 	vMultScal(1/sqrt(tmp),prop_0,prop_0);
@@ -1894,6 +1954,9 @@ void InitVariables(void)
 	orient_used=false;
 	directory="";
 	lambda=TWO_PI;
+	beam_center_used=false;
+	deprecated_bc_used=false;
+	vInit(beam_center_0);
 	// initialize ref_index of scatterer
 	Nmat=Nmat_given=1;
 	ref_index[0]=1.5;
@@ -2075,6 +2138,9 @@ void VariablesInterconnect(void)
 	/* TO ADD NEW INTERACTION FORMULATION
 	 * If the new Green's tensor is non-symmetric (which is very unlikely) add it to the definition of reduced_FFT
 	 */
+	
+	if (deprecated_bc_used && beam_center_used) LogError(ONE_POS,"Beam center coordinates can not be specified as "
+		"arguments to both '-beam' and '-beam_center'. Use only the latter.");
 	if (calc_Csca || calc_vec) all_dir = true;
 	// by default, one of the scattering options is activated
 	if (store_scat_grid || phi_integr) scat_grid = true;
@@ -2430,6 +2496,7 @@ void PrintInfo(void)
 		fprintf(logfile,"Volume-equivalent size parameter: "GFORM"\n",ka_eq);
 		// log incident beam and polarization
 		fprintf(logfile,"\n---In laboratory reference frame:---\nIncident beam: %s\n",beam_descr);
+		fprintf(logfile,"Incident beam center position: "GFORMDEF3V"\n",COMP3V(beam_center_0));
 		fprintf(logfile,"Incident propagation vector: "GFORMDEF3V"\n",COMP3V(prop_0));
 		if (beamtype==B_DIPOLE) fprintf(logfile,"(dipole orientation)\n");
 		else { // polarizations are not shown for dipole incident field
@@ -2454,8 +2521,7 @@ void PrintInfo(void)
 			if (alph_deg!=0 || bet_deg!=0 || gam_deg!=0) {
 				fprintf(logfile,"Particle orientation (deg): alpha="GFORMDEF", beta="GFORMDEF", gamma="GFORMDEF"\n\n"
 					"---In particle reference frame:---\n",alph_deg,bet_deg,gam_deg);
-				if (beam_asym) fprintf(logfile,"Incident Beam center position: "GFORMDEF3V"\n",
-					beam_center[0],beam_center[1],beam_center[2]);
+				fprintf(logfile,"Incident beam center position: "GFORMDEF3V"\n",COMP3V(beam_center_0));
 				fprintf(logfile,"Incident propagation vector: "GFORMDEF3V"\n",COMP3V(prop));
 				if (beamtype==B_DIPOLE) fprintf(logfile,"(dipole orientation)\n");
 				else { // polarizations are not shown for dipole incident field
diff --git a/src/vars.c b/src/vars.c
index 7d54d7bb..347772b9 100644
--- a/src/vars.c
+++ b/src/vars.c
@@ -63,12 +63,16 @@ double propAlongZ;  // equal 0 for general incidence, and +-1 for incidence alon
 bool rectDip;       // whether using rectangular-cuboid (non-cubical) dipoles
 
 // 3D vectors (in particle reference frame)
-double prop_0[3],prop[3];     // incident direction (in laboratory and particle reference frame)
-double incPolX[3],incPolY[3]; // incident polarizations (in particle RF)
+double prop[3];               // incident direction
+double incPolX[3],incPolY[3]; // incident polarizations
 double beam_center[3];        // coordinates of the beam center
 double box_origin_unif[3];    /* coordinates of the center of the first dipole in the local computational box (after
                                  uniform distribution of non-void dipoles among all processors) */
 
+// 3D vectors (in laboratory reference frame)
+double prop_0[3];             // incident direction 
+double beam_center_0[3];      // coordinates of the beam center
+
 // file info
 const char * restrict directory; // directory to save data in
 FILE * restrict logfile;         // file where all the information about the run is saved
diff --git a/src/vars.h b/src/vars.h
index a7fcc267..356f1c7b 100644
--- a/src/vars.h
+++ b/src/vars.h
@@ -47,7 +47,7 @@ extern bool prognosis,yzplane,scat_plane,store_mueller,all_dir,scat_grid,phi_int
 extern double propAlongZ;
 
 // 3D vectors
-extern double prop_0[3],prop[3],incPolX[3],incPolY[3],beam_center[3],box_origin_unif[3];
+extern double prop_0[3],prop[3],incPolX[3],incPolY[3],beam_center_0[3],beam_center[3],box_origin_unif[3];
 
 // file info
 extern const char * restrict directory;