| 99 | |
| 100 | |
| 101 | def write_output(c): |
| 102 | e = np.zeros((n + 1, m + 1)) |
| 103 | k = 0 |
| 104 | edata = '' |
| 105 | for i in range(2, n + 1): |
| 106 | erow = '' |
| 107 | for j in range(0, m + 1): |
| 108 | e[i][j] = c[k] |
| 109 | erow += 'e[{i:d}][{j:d}] = {val:0.16g}; '.format(val=e[i][j], i=i, j=j) |
| 110 | k += 1 |
| 111 | edata += erow + '\n' |
| 112 | f1 = c[k] |
| 113 | g1 = c[k + 1] |
| 114 | g2 = c[k + 2] |
| 115 | |
| 116 | template = textwrap.dedent( |
| 117 | """ |
| 118 | double {name:s}Class::viscosity_Trho(double T, double rho) |
| 119 | {{ |
| 120 | // This function was generated by fitting REFPROP ECS data |
| 121 | // to the functional form of Vogel, 1998 (propane viscosity) |
| 122 | // The script entitled dev/fit_avoid_ECS.py was used to make this |
| 123 | // function. The mean absolute error of the fit is equal to |
| 124 | // {MAE:g} % |
| 125 | |
| 126 | double delta_0, sum, DELTA_H_eta, e_k, sigma, tau, delta; |
| 127 | double e[{n:d}+1][{m:d}+1]; |
| 128 | |
| 129 | tau = T/reduce.T; //[Opposite to normal definition] |
| 130 | delta = rho/reduce.rho; |
| 131 | |
| 132 | // Load the coefficients |
| 133 | double f1 = {f1:0.16g}, g1 = {g1:0.16g}, g2 = {g2:0.16g}; |
| 134 | for (int i=0;i<={n:d};i++){{ for(int j=0;j<={m:d};j++){{ e[i][j]=0.0; }} }} |
| 135 | {edata:s} |
| 136 | delta_0=g1*(1+g2*sqrt(tau)); //[no units] |
| 137 | sum=0; |
| 138 | for (int i=2;i<={n:d};i++){{ |
| 139 | for (int j=0;j<={m:d};j++){{ |
| 140 | sum += e[i][j]*pow(delta,i)/pow(tau,j); |
| 141 | }} |
| 142 | }} |
| 143 | DELTA_H_eta = sum + f1*(delta/(delta_0-delta)-delta/delta_0); //[Pa-s] |
| 144 | |
| 145 | try{{ |
| 146 | // Get the ECS params for the fluid if it has them |
| 147 | ECSParams(&e_k,&sigma); |
| 148 | }} |
| 149 | catch(NotImplementedError) |
| 150 | {{ |
| 151 | throw ValueError(format("Your fluid does not implement ECSParams")); |
| 152 | }} |
| 153 | |
| 154 | return viscosity_dilute(T,e_k,sigma) + DELTA_H_eta; |
| 155 | }} |
| 156 | """ |
| 157 | ) |
| 158 | |