0001 void fluctuation(mode,krho,m,n,Y,krho_s1,krho_s2,krho_s3,krho_s4,krho_s5,krho_s6, m_s1,m_s2,m_s3,m_s4,m_s5, n_s1,n_s2,n_s3,n_s4,n_s5,raymagnetic,rayprofiles,param,equilparam,Npar,ap,cnorm,ran2,res,reset)
0002 int mode;
0003 double *krho;
0004 double *m;
0005 double *n;
0006 double Y[];
0007 double krho_s1,krho_s2,krho_s3,krho_s4,krho_s5,krho_s6;
0008 double m_s1,m_s2,m_s3,m_s4,m_s5;
0009 double n_s1,n_s2,n_s3,n_s4,n_s5;
0010 raymagnetic_format raymagnetic;
0011 rayprofiles_format rayprofiles;
0012 rayparam_format param;
0013 equilparam_format equilparam;
0014 double Npar;
0015 double ap,cnorm;
0016 double clum=CLUM;
0017 double ran2;
0018 long int *reset;
0019 double *res;
0020 double pi = 3.141592653526;
0021 {
0022 static double ran2(idum)
0023 long int *idum;
0024 {
0025
0026 long int j;
0027 long int k;
0028 static long int idum2 = 123456789;
0029 static long int iy = 0;
0030 static long int iv[32];
0031 double temp;
0032
0033 double eps = 1.2e-7,rnmx,am,ndiv;
0034 long int im1 = 2147483563,im2 = 2147483399,imm1,ia1 = 40014,ia2 = 40692,iq1 = 53668,iq2 = 52774;
0035 long int ir1 = 12211,ir2 = 3791,ntab = 32;
0036
0037 imm1 = im1 - 1L;
0038 am = (double)1.0/(double)im1;
0039 ndiv = (double)1.0+(double)imm1/(double)ntab;
0040 rnmx = (double)1.0 - eps;
0041
0042 if (*idum <=0L) {
0043 if(-(*idum) < 1L) *idum = 1L;
0044 else *idum = -(*idum);
0045 idum2 = (*idum);
0046 for (j = ntab+7L;j>=0L;j--) {
0047 k = (long int)(*idum)/iq1;
0048 *idum = ia1*(*idum-k*iq1)-k*ir1;
0049 if(*idum<0L) *idum += im1;
0050 if (j < ntab) iv[j] = *idum;
0051 }
0052 iy = iv[0];
0053 }
0054
0055 k = (long int)(*idum)/iq1;
0056 *idum = ia1*(*idum-k*iq1)-k*ir1;
0057 if (*idum < 0L) *idum += im1;
0058 k = (long int)(idum2/iq2);
0059 idum2 = ia2*(idum2-k*iq2)-k*ir2;
0060 if (idum2 < 0L) idum2 += im2;
0061 j = (long)(iy/(long int)ndiv);
0062 iy = iv[j]-idum2;
0063 iv[j] = *idum;
0064 if(iy < 1) iy += imm1;
0065 if ((temp = am*(double)iy) > rnmx) return(rnmx);
0066 else return(temp);
0067 }
0068 srand(time(NULL));
0069 reset=rand();
0070 res[0] = (double)2.0*pi*ran2(reset);
0071 ap = equilparam.ap[0];
0072 cnorm = clum/ap/Y[6];
0073
0074 krho_s1=raymagnetic.Bsn[0]*(Npar/cnorm)*(raymagnetic.salpha[0]/raymagnetic.calpha[0]);
0075 krho_s2=(raymagnetic.BPHIn[0]/raymagnetic.calpha[0])*(Y[4]/raymagnetic.rn[0]);
0076 krho_s3=-raymagnetic.BPHIn[0]*Y[3]*raymagnetic.gradrhon[0]*(raymagnetic.salpha[0]/raymagnetic.calpha[0]);
0077 krho_s4=-raymagnetic.Bsn[0]*Y[5]/raymagnetic.Rn[0];
0078 krho_s5=raymagnetic.Bsn[0]*(Npar/cnorm)*(raymagnetic.salpha[0]/raymagnetic.calpha[0]);
0079 krho_s6=-Y[3]*raymagnetic.gradrhon[0];
0080
0081 krho[0]=(krho_s1+((krho_s2+krho_s3+krho_s4)*sin(res[0]))-((krho_s5+krho_s6)*cos(res[0])))/raymagnetic.gradrhon[0];
0082
0083 m_s1=(raymagnetic.Bsn[0]/raymagnetic.calpha[0])*(Npar/cnorm);
0084 m_s2=raymagnetic.BPHIn[0]*(Y[4]/raymagnetic.rn[0])*(raymagnetic.salpha[0]/raymagnetic.calpha[0]);
0085 m_s3=-(raymagnetic.BPHIn[0]/raymagnetic.calpha[0])*Y[3]*raymagnetic.gradrhon[0];
0086 m_s4=(raymagnetic.Bsn[0]/raymagnetic.calpha[0])*(Npar/cnorm);
0087 m_s5=-Y[4]/raymagnetic.rn[0];
0088
0089 m[0]=(m_s1+((m_s2+m_s3)*sin(res[0]))-((m_s4+m_s5)*cos(res[0])))*raymagnetic.rn[0];
0090
0091 n_s1=raymagnetic.BPHIn[0]*(Npar/cnorm);
0092 n_s2=raymagnetic.Bsn[0]*Y[3]*raymagnetic.gradrhon[0];
0093 n_s3=-raymagnetic.Bsn[0]*(Y[4]/raymagnetic.rn[0])*raymagnetic.salpha[0];
0094 n_s4=raymagnetic.BPHIn[0]*(Npar/cnorm);
0095 n_s5=-Y[5]/raymagnetic.Rn[0];
0096
0097 n[0]=(n_s1+((n_s2+n_s3)*sin(res[0]))-((n_s4+n_s5)*cos(res[0])))*raymagnetic.Rn[0];
0098 }