firstordercollop_dke_ft_c

C/C++ source

0001 #include <stdio.h>
0002 #include <stdlib.h>
0003 #include <math.h>
0004 #include "Fokker_Planck_DKE_loop.h"
0005 
0006 
0007 int firstordercollop_dke_ft_c2(struct Tensor f1,
0008 struct Tensor betath_ref,
0009 struct Tensor masku,struct Tensor maskl, 
0010 struct Tensor dpn, 
0011 struct Tensor xpn, struct Tensor xpn2, struct Tensor xpn3, struct Tensor xsigma, struct Tensor xgamma,
0012 struct Tensor xz, struct Tensor xz2,
0013 struct Tensor pn, struct Tensor pn2, struct Tensor pn3,  
0014 struct Tensor sigma, struct Tensor gamma,
0015 struct Tensor z, struct Tensor z2, 
0016 struct Tensor xJ1, struct Tensor xJ2, struct Tensor xJ3, struct Tensor xJ4,
0017 struct Tensor J1, struct Tensor J2, struct Tensor J3, struct Tensor J4,struct
0018 Tensor xTe_norm,struct Tensor *I1,struct Tensor *I2)
0019 
0020 {
0021 /* printf("FIRSTORDERCOLLOP\n");*/
0022 int i,j,sdim[2];
0023 struct Tensor A;
0024 struct Tensor B1,B2,B3,B5; /*A,B are intermediate tensor variables*/
0025 struct Tensor xf1;
0026 
0027 
0028 /*Creating xf1*/
0029 
0030 if(f1.ndims !=2 && (f1.sizedim[0] !=1 || f1.sizedim[1]!=1)){fprintf(stdout,"error:f1 should be a vector");return -1;}
0031 else if (f1.ndims ==2 && f1.sizedim[0] ==1){i=1;}
0032 else if (f1.ndims ==2 && f1.sizedim[1] ==1){i=0;}
0033 
0034 
0035 sdim[0]=f1.sizedim[i];
0036 sdim[1]=f1.sizedim[i];
0037 Tensor_Create(2,sdim,&xf1);
0038 for(i=0;ifor(j=0;j/*ith line of xf1 == f1.vals[i]*/
0041  }
0042 }
0043 
0044 /*               */
0045 
0046 /*Computing tensors products for I1*/
0047 
0048 /*1st line*/
0049 
0050 product(xf1,masku,&A);
0051 product(xpn3,A,&A);/*A=xf1;*masku.*xpn3*/
0052 division(A,xgamma,&B2);
0053 integral_dke_jd(dpn,B2,1,&(*I1));
0054 Scalar_Tensor_product(1/3.0,(*I1),&(*I1));
0055 division((*I1),xTe_norm,&(*I1));/*Could rewrite xTe_norm as scalar and use
0056 Scalar_Tensor_product instead*/
0057 
0058 /*2nd line*/
0059 
0060 integral_dke_jd(dpn,A,1,&B2);
0061 product(gamma,B2,&B2);
0062 Scalar_Tensor_product(2.0/3.0,B2,&B2);
0063 division(B2,xTe_norm,&B2);
0064 Tensor_Sum((*I1),B2,-1.0,&(*I1));
0065 
0066 /*3rd line*/
0067 
0068 product(A,xpn2,&B2);
0069 division(B2,xgamma,&B2);
0070 integral_dke_jd(dpn,B2,1,&B3);
0071 product(B3,gamma,&B3);
0072 Scalar_Tensor_product(1.0/5.0,B3,&B3);
0073 division(B3,xTe_norm,&B3);
0074 division(B3,xTe_norm,&B3);
0075 Tensor_Sum((*I1),B3,1.0,&(*I1));
0076 
0077 /*5th line*/
0078 
0079 product(A,xJ2,&B2);
0080 division(B2,xgamma,&B2);
0081 division(B2,xz2,&B2);
0082 integral_dke_jd(dpn,B2,1,&B3);
0083 product(B3,gamma,&B3);
0084 division(B3,xTe_norm,&B3);
0085 Tensor_Sum((*I1),B3,-1.0,&(*I1));
0086 
0087 /*6th line*/
0088 
0089 struct Tensor B4;
0090 
0091 Tensor_Create(xz2.ndims,xz2.sizedim, &B4);
0092 for(i=0;ifor(i=0;ifor(i=0;i/*7th line*/
0110 /* fprintf(stdout,"line 7=======================\n");*/
0111 product(A,xJ3,&B3);
0112 division(B3,xz,&B3);
0113 division(B3,xgamma,&B3);
0114 integral_dke_jd(dpn,B3,1,&B2);
0115 product(B2,gamma,&B2);
0116 division(B2,betath_ref,&B2);
0117 division(B2,xTe_norm,&B2);
0118 division(B2,betath_ref,&B2);
0119 division(B2,xTe_norm,&B2);
0120 for(i=0;i/*8th line*/
0124 /* fprintf(stdout,"line 8=======================\n");*/
0125 product(A,xJ1,&B3);
0126 division(B3,xz2,&B3);
0127 division(B3,xgamma,&B3);
0128 integral_dke_jd(dpn,B3,1,&B2);
0129 product(B2,gamma,&B2);
0130 division(B2,xTe_norm,&B2);
0131 for(i=0;i/*10th line*/
0135 /* fprintf(stdout,"line 10=======================\n");*/
0136 product(A,xJ4,&B3);
0137 division(B3,xz,&B3);
0138 division(B3,xgamma,&B3);
0139 integral_dke_jd(dpn,B3,1,&B2);
0140 product(B2,gamma,&B2);
0141 division(B2,xTe_norm,&B2);
0142 division(B2,betath_ref,&B2);
0143 division(B2,xTe_norm,&B2);
0144 division(B2,betath_ref,&B2);
0145 for(i=0;i/*9th line*/
0149 
0150 product(xf1,masku,&B3);
0151 product(B3,xpn,&B3);
0152 product(xgamma,xsigma,&B4);
0153 division(B4,xz,&B4);
0154 for(i=0;i/*4th line*/
0163 
0164 division(xsigma,xz,&A);
0165 Tensor_Sum(xgamma,A,-1.0,&A);
0166 product(A,xpn,&A);
0167 division(A,xgamma,&A);
0168 product(A,masku,&A);
0169 product(A,xf1,&A);
0170 integral_dke_jd(dpn,A,1,&B2);
0171 product(gamma,pn2,&B4);
0172 
0173 Tensor_Sum((*I1),B2,1.0,&(*I1));
0174 
0175 /*=========================================================================*/
0176 /*Computing tensors products for I2*/
0177 
0178 product(xf1,maskl,&A);
0179 
0180 /*1st line*/
0181 
0182 division(A,xgamma,&B3);
0183 integral_dke_jd(dpn,B3,1,&B2);
0184 Tensor_Create(B2.ndims, B2.sizedim, &(*I2));
0185 Tensor_Fill(I2,0.0);
0186 division(B2,xTe_norm,&B2);
0187 for(i=0;i/*2nd line*/
0191 
0192 integral_dke_jd(dpn,A,1,&B2);
0193 division(gamma,xTe_norm,&B3);
0194 for(i=0;ifor(i=0;i/*3rd line*/
0204 
0205 division(A,xgamma,&B4);
0206 integral_dke_jd(dpn,B4,1,&B3);
0207 division(sigma,z,&B2);
0208 Tensor_Sum(gamma,B2,-1.0,&B2);
0209 division(B2,pn2,&B2);
0210 product(B3,B2,&B2);
0211 Tensor_Sum((*I2),B2,1.0,&(*I2));
0212 
0213 
0214 /*4th line*/
0215 
0216 integral_dke_jd(dpn,A,1,&B2);
0217 product(J2,B2,&B2);
0218 division(B2,z2,&B2);
0219 division(B2,xTe_norm,&B2);
0220 Tensor_Sum((*I2),B2,-1.0,&(*I2));
0221 
0222 
0223 /*5th line*/
0224 
0225 product(xgamma,xpn2,&B4);
0226 for(i=0;ifor(i=0;ifor(i=0;i/*6th line*/
0242 
0243 division(J3,z,&B4);
0244 division(B4,betath_ref,&B4);division(B4,betath_ref,&B4);
0245 division(B4,xTe_norm,&B4);division(B4,xTe_norm,&B4);
0246 for(i=0;ifor(i=0;ifor(i=0;i/*smth wrong with the sign of B2;changed to -1*/
0259 
0260 /*7th line*/
0261 
0262 product(A,xpn2,&B3);
0263 division(B3,xgamma,&B3);
0264 integral_dke_jd(dpn,B3,1,&B2);
0265 product(gamma,sigma,&B4);
0266 division(B4,z,&B4);
0267 for(i=0;ireturn 0;}