//function defines the equations of GR polarized radiative transfer in polarized coordinates
int trans (doub llog, const doub yyy[], doub ff[], void *pas) {
	//sinh transformation of affine parameter t allows to convert scales close to the BH and far from the BH to a single uniform grid
     doub det=-3.*cosh(llog),
          t=1.+ 3.* sinh(llog)/r0;
	 //evaluation of plasma emissivities/absorptivities/rotativities at a given point
#include "evalpointzero.cpp"

    //local auxiliary variables
    doub cosqs=cos(yyy[2]), 
		 sinqs=sin(yyy[2]), 
		 frsq=fr*fr, 
		 frcu=frsq*fr, 
		 sinths=sin(yyy[3]), 
		 cosths=cos(yyy[3]);

	//radiative transfer equations in polarized coordinates
    //this piece of code is autogenerated in Mathematica
    ff[0]=(det*(jIc - frcu*(aIc*yyy[0] + (aVc*cosths +aQc*cos2k*cosqs*sinths + aQc*sin2k*sinqs*sinths)*yyy[1])))/frsq;
    ff[2]=(det*(cosqs*sin2k - cos2k*sinqs)*(jQc - aQc*frcu*yyy[0]) - det*frcu*(cos2k*cosqs*cosths*rQc + cosths*rQc*sin2k*sinqs -rVc*sinths)*yyy[1])/(frsq*sinths*yyy[1]);
    ff[1]= (det*((cos2k*cosqs + sin2k*sinqs)*sinths*(jQc - aQc*frcu*yyy[0]) + cosths*(jVc -aVc*frcu*yyy[0]) - aIc*frcu*yyy[1]))/frsq;
    ff[3]=(det*(-(jVc*sinths) + aVc*frcu*sinths*yyy[0] + cos2k*cosqs*cosths*(jQc - aQc*frcu*yyy[0]) + cosths*sin2k*sinqs*(jQc 
	    - aQc*frcu*yyy[0]) + cosqs*frcu*rQc*sin2k*yyy[1] -cos2k*frcu*rQc*sinqs*yyy[1]))/(frsq*yyy[1]);
    ff[4]= det*(-fr*xaIc*yyy[4] + xjIc/frsq);

    //Check for NANs
    if(ff[0]!=ff[0]) {
    	printf("Error in radiative transfer \n Exiting \n"); 
	    exit(1);
    };
    return(0);
}