HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
utoprim_jon_eos.c
Go to the documentation of this file.
1 
13 /*
14  This file must contain the equation of state
15  in two different forms:
16 
17  p(rho0,w)
18  p(rho0,u)
19 
20  and also the derivatives
21 
22  dp/du
23  dp/drho0
24 
25  which are required only to evaluate the
26  matrix dU/dprim in Utoprim_5D; these
27  are not required for Utoprim_1D.
28 
29  current equation of state is a GAMMA-law gas.
30 
31 */
32 
33 
34 
35 
36 
37 // 1/ (dwmrho0/dp)
38 // used for dP/dW and dP/dvsq
39 //static FTYPE compute_idwmrho0dp_old(FTYPE rho0, FTYPE wmrho0)
40 //{
41 // static FTYPE dpressuredu_rho0_u(FTYPE rho0, FTYPE u);
42 // static FTYPE u_wmrho0(FTYPE rho0, FTYPE wmrho0);
43 // FTYPE u,dpdu;
44 // u=u_wmrho0(rho0, wmrho0);
45 // dpdu=dpressuredu_rho0_u(rho0, u);
46 //
47 //return(1.0/(1.0/dpdu+1.0));
48 //}
49 
50 
51 // static declarations
52 static FTYPE pressure_Wp_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq);
53 static FTYPE pressure_W_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) ;
54 static FTYPE wmrho0_compute_vsq(FTYPE Wp, FTYPE D, FTYPE vsq, FTYPE gamma,FTYPE gammasq);
55 static FTYPE pressure_W_vsq_scn(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) ;
56 static FTYPE pressure_W_vsq_old(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) ;
57 static FTYPE wmrho0_compute_utsq(FTYPE Wp, FTYPE D, FTYPE utsq, FTYPE gamma,FTYPE gammasq);
58 static FTYPE dpdWp_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq);
59 static FTYPE dSsdWp_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq);
60 static FTYPE wmrho0_compute_vsq(FTYPE Wp, FTYPE D, FTYPE vsq, FTYPE gamma,FTYPE gammasq);
61 static FTYPE wmrho0_compute_utsq(FTYPE Wp, FTYPE D, FTYPE utsq, FTYPE gamma,FTYPE gammasq);
62 static FTYPE dpdW_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq);
63 static FTYPE dpdvsq_calc_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq);
64 static FTYPE dSsdW_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq);
65 static FTYPE dSsdvsq_calc_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq);
66 static FTYPE dpdvsq_calc_scn(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq);
67 static FTYPE wmrho0_compute_vsq(FTYPE Wp, FTYPE D, FTYPE vsq, FTYPE gamma,FTYPE gammasq);
68 static FTYPE wmrho0_compute_utsq(FTYPE Wp, FTYPE D, FTYPE utsq, FTYPE gamma,FTYPE gammasq);
69 static FTYPE wmrho0_compute_vsq(FTYPE Wp, FTYPE D, FTYPE vsq, FTYPE gamma,FTYPE gammasq);
70 
71 
72 /*
73 
74  pressure as a function of W, vsq, and D:
75 
76 
77 */
78 
79 //#define GAMMA (gam)
80 
81 #define OLDCALCJON 0
82 
83 //#define CRAZYDEBUG 1
84 
85 #if(OLDCALCJON)
86 static FTYPE pressure_W_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
87 {
88  FTYPE gtmp;
89  FTYPE answer;
90 
91  gtmp = 1. - vsq;
92 
93  answer=(GAMMA - 1.) * ( W * gtmp - D * sqrt(gtmp) ) / GAMMA ;
94 
95  return(answer);
96 
97 }
98 #else
99 static FTYPE pressure_W_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
100 {
101 
102  return(pressure_Wp_vsq(whicheos,EOSextra, W-D,D,vsq));
103 
104 }
105 #endif
106 
107 
108 
109 #if(OLDCALCJON)
110 static FTYPE pressure_Wp_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
111 {
112 
113 
114  return(pressure_W_vsq(whicheos,EOSextra, Wp+D, D, vsq));
115 
116 }
117 
118 #else
119 
120 
121 static FTYPE pressure_Wp_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
122 {
123  FTYPE gtmp;
124  FTYPE gamma,gammasq;
125  FTYPE rho0,wmrho0;
126  FTYPE pressure_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
127 
128 
129 
130  gamma=1.0/sqrt(1.0-vsq);
131  gammasq=gamma*gamma;
132  wmrho0=wmrho0_compute_vsq(Wp, D, vsq, gamma,gammasq);
133  rho0=D/gamma;
134 
135  // wmrho0=u+p and for ideal gas p=(GAMMA-1) u, so p (GAMMA/(GAMMA-1)) = wmrho0
136  // return( wmrho0*IGAMMAR );
137  return(pressure_wmrho0(whicheos,EOSextra,rho0,wmrho0));
138 
139 
140  // return(pressure_W_vsq_scn(Wp+D,D,vsq));
141 
142  // return(pressure_W_vsq_old(Wp+D,D,vsq));
143 }
144 
145 static FTYPE pressure_Wp_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq)
146 {
147  FTYPE gtmp;
148  FTYPE gamma,gammasq;
149  FTYPE rho0,wmrho0;
150  FTYPE pressure_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
151 
152 
153  gammasq=1.0+utsq;
154  gamma=sqrt(gammasq);
155  wmrho0=wmrho0_compute_utsq(Wp, D, utsq, gamma,gammasq);
156  rho0=D/gamma;
157 
158  // wmrho0=u+p and for ideal gas p=(GAMMA-1) u, so p (GAMMA/(GAMMA-1)) = wmrho0
159  return(pressure_wmrho0(whicheos,EOSextra,rho0,wmrho0));
160 
161 
162 }
163 
164 static FTYPE Ss_Wp_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq)
165 {
166  FTYPE gtmp;
167  FTYPE gamma,gammasq;
168  FTYPE rho0,wmrho0;
169  FTYPE compute_specificentropy_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
170 
171 
172  gammasq=1.0+utsq;
173  gamma=sqrt(gammasq);
174  wmrho0=wmrho0_compute_utsq(Wp, D, utsq, gamma,gammasq);
175  rho0=D/gamma;
176 
177  // wmrho0=u+p and for ideal gas p=(GAMMA-1) u, so p (GAMMA/(GAMMA-1)) = wmrho0
178  return(compute_specificentropy_wmrho0(whicheos,EOSextra,rho0,wmrho0));
179 
180 
181 }
182 
183 
184 
185 #endif
186 
187 
188 static FTYPE pressure_W_vsq_old(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
189 {
190  FTYPE gtmp;
191 
192 
193  gtmp = 1. - vsq;
194 
195  return( (gam - 1.) * ( W * gtmp - D * sqrt(gtmp) ) / gam );
196 
197 }
198 
199 
200 static FTYPE pressure_W_vsq_scn(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
201 {
202  FTYPE gtmp;
203  FTYPE gamma;
204  FTYPE wmrho0;
205  FTYPE rho0;
206  FTYPE pressure_wmrho0(int whicheos,FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
207 
208 
209  gamma=1.0/sqrt(1.0-vsq);
210  rho0=D/gamma;
211 
212  gtmp = 1. - vsq;
213  wmrho0=W * gtmp - D * sqrt(gtmp);
214 
215  return(pressure_wmrho0(whicheos, EOSextra, rho0, wmrho0));
216 
217 
218 }
219 
220 
221 
222 /*
223 
224  partial derivative of pressure with respect to W holding vsq fixed
225 
226  // dp/dW|v^2 = dp/drho|wmrho0 drho/dW|v^2 + dp/dwmrho0|rho dwmrho0/dW|v^2
227 
228  */
229 #if(OLDCALCJON)
230 static FTYPE dpdW_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
231 {
232 
233  return( (GAMMA - 1.) * (1. - vsq) / GAMMA ) ;
234 
235 }
236 
237 static FTYPE dpdWp_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
238 {
239  FTYPE W=Wp+D;
240 
241  return( (GAMMA - 1.) * (1. - vsq) / GAMMA ) ;
242 
243 }
244 
245 #else
246 
247 static FTYPE dpdW_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
248 {
249 
250 
251  return(dpdWp_calc_vsq(whicheos,EOSextra, W-D, D, vsq));
252 
253 }
254 
255 
258 static FTYPE dpdWp_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
259 {
260  FTYPE compute_idwmrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
261  FTYPE compute_idrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
262  FTYPE rho0,wmrho0;
263  FTYPE idwmrho0dp;
264  FTYPE gamma,gammasq;
265  FTYPE dwmrho0dW;
266  FTYPE dpdW;
267  FTYPE drho0dW,idrho0dp;
268 
269  gamma=1.0/sqrt(1.0-vsq);
270  gammasq=gamma*gamma;
271  wmrho0=wmrho0_compute_vsq(Wp, D, vsq, gamma, gammasq);
272 
273  rho0=D/gamma;
274  // holding v^2 fixed
275  idwmrho0dp=compute_idwmrho0dp(whicheos, EOSextra, rho0, wmrho0);
276  dwmrho0dW = 1.0-vsq;
277 
278  drho0dW = 0.0; // because \rho=D/\gamma
279  idrho0dp = compute_idrho0dp(whicheos, EOSextra, rho0, wmrho0); // so not really needed
280 
281  dpdW = drho0dW *idrho0dp + dwmrho0dW *idwmrho0dp;
282 
283  return( dpdW ) ;
284 
285 }
286 
289 static FTYPE dpdWp_calc_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq)
290 {
291  FTYPE compute_idwmrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
292  FTYPE compute_idrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
293  FTYPE rho0,wmrho0;
294  FTYPE idwmrho0dp;
295  FTYPE gamma,gammasq;
296  FTYPE dwmrho0dW;
297  FTYPE dpdW;
298  FTYPE drho0dW,idrho0dp;
299 
300  gammasq=1.0+utsq;
301  gamma=sqrt(gammasq);
302  wmrho0=wmrho0_compute_utsq(Wp, D, utsq, gamma, gammasq);
303 
304  rho0=D/gamma;
305  // holding utsq fixed
306  idwmrho0dp=compute_idwmrho0dp(whicheos,EOSextra, rho0, wmrho0);
307  dwmrho0dW = 1.0/gammasq;
308 
309  drho0dW = 0.0; // because \rho=D/\gamma
310  idrho0dp = compute_idrho0dp(whicheos,EOSextra, rho0, wmrho0); // so not really needed
311 
312  dpdW = drho0dW *idrho0dp + dwmrho0dW *idwmrho0dp;
313 
314  return( dpdW ) ;
315 
316 }
317 
318 
319 
320 
321 
324 static FTYPE dSsdWp_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
325 {
326  FTYPE compute_dspecificSdrho_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
327  FTYPE compute_dspecificSdwmrho0_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
328  FTYPE rho0,wmrho0;
329  FTYPE idwmrho0dSs;
330  FTYPE gamma,gammasq;
331  FTYPE dwmrho0dW;
332  FTYPE dSsdW;
333  FTYPE drho0dW,idrho0dSs;
334 
335  gamma=1.0/sqrt(1.0-vsq);
336  gammasq=gamma*gamma;
337  wmrho0=wmrho0_compute_vsq(Wp, D, vsq, gamma, gammasq);
338 
339  rho0=D/gamma;
340  // holding v^2 fixed
341  idwmrho0dSs=compute_dspecificSdwmrho0_wmrho0(whicheos, EOSextra, rho0, wmrho0);
342  dwmrho0dW = 1.0-vsq;
343 
344  drho0dW = 0.0; // because \rho=D/\gamma
345  idrho0dSs = compute_dspecificSdrho_wmrho0(whicheos, EOSextra, rho0, wmrho0); // so not really needed
346 
347  dSsdW = drho0dW *idrho0dSs + dwmrho0dW *idwmrho0dSs;
348 
349  return( dSsdW ) ;
350 
351 }
352 
355 static FTYPE dSsdWp_calc_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq)
356 {
357  FTYPE compute_dspecificSdrho_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
358  FTYPE compute_dspecificSdwmrho0_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
359  FTYPE rho0,wmrho0;
360  FTYPE idwmrho0dSs;
361  FTYPE gamma,gammasq;
362  FTYPE dwmrho0dW;
363  FTYPE dSsdW;
364  FTYPE drho0dW,idrho0dSs;
365 
366  gammasq=1.0+utsq;
367  gamma=sqrt(gammasq);
368  wmrho0=wmrho0_compute_utsq(Wp, D, utsq, gamma, gammasq);
369 
370  rho0=D/gamma;
371  // holding utsq fixed
372  idwmrho0dSs=compute_dspecificSdwmrho0_wmrho0(whicheos,EOSextra, rho0, wmrho0);
373  dwmrho0dW = 1.0/gammasq;
374 
375  drho0dW = 0.0; // because \rho=D/\gamma
376  idrho0dSs = compute_dspecificSdrho_wmrho0(whicheos,EOSextra, rho0, wmrho0); // so not really needed
377 
378  dSsdW = drho0dW *idrho0dSs + dwmrho0dW *idwmrho0dSs;
379 
380  return( dSsdW ) ;
381 
382 }
383 
384 
385 
386 
387 
388 
389 
390 #endif
391 
392 
393 
394 
395 
396 
397 static FTYPE dpdW_calc_vsq_scn(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
398 {
399 
400  return(dpdW_calc_vsq(whicheos,EOSextra,W,D,vsq));
401 
402 }
403 
404 /*
405 
406  partial derivative of pressure with respect to vsq holding W fixed
407 
408 
409 */
410 #if(OLDCALCJON)
411 static FTYPE dpdvsq_calc(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
412 {
413  FTYPE outval;
414  outval = (GAMMA - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / GAMMA ;
415 
416 
417  return(outval);
418 }
419 
420 #else
421 static FTYPE dpdvsq_calc(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
422 {
423 
424  // old version
425  // return( (GAMMA - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / GAMMA );
426 
427  return(dpdvsq_calc_Wp(whicheos,EOSextra, W-D, D, vsq));
428 }
429 #endif
430 
431 // GODMARK: probably should write "clean" version of this, but maybe doesn't matter
432 
433 
434 #if(OLDCALCJON)
435 static FTYPE dpdvsq_calc_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
436 {
437 
438  return(dpdvsq_calc_scn(whicheos,EOSextra, Wp+D, D, vsq));
439 }
440 
441 #else
442 
443 
445 static FTYPE dpdvsq_calc_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
446 {
447  FTYPE outval;
448  FTYPE gamma,W;
449  FTYPE compute_idwmrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
450  FTYPE compute_idrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
451  FTYPE rho0,wmrho0;
452  FTYPE dwmrho0dvsq,idwmrho0dp;
453  FTYPE drho0dvsq,idrho0dp;
454  FTYPE gammasq;
455  FTYPE dpdvsq;
456 
457 
458  gamma=1.0/sqrt(1.0-vsq);
459  gammasq=gamma*gamma;
460  wmrho0=wmrho0_compute_vsq(Wp, D, vsq, gamma, gammasq);
461 
462  rho0=D/gamma;
463 
464  // return( (gam - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / gam ) ;
465  // outval = (gam - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / gam ;
466 
467  W = Wp+D;
468 
469  // outval = (D*gamma*0.5-W)*idwmrho0dp; // old code
470 
471  idwmrho0dp=compute_idwmrho0dp(whicheos,EOSextra, rho0, wmrho0);
472  dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp);
473 
474  drho0dvsq = -D*gamma*0.5; // because \rho=D/\gamma
475  idrho0dp = compute_idrho0dp(whicheos,EOSextra, rho0, wmrho0);
476 
477  dpdvsq = drho0dvsq *idrho0dp + dwmrho0dvsq *idwmrho0dp;
478 
479  // if( outval > 0. ) {
480  // dualfprintf(fail_file,"outval = %21.15g , D = %21.15g , vsq = %21.15g, W = %21.15g \n",
481  // outval, D, vsq, W );
482  // }
483 
484  return(dpdvsq);
485 }
486 
487 
489 static FTYPE dpdvsq_calc2_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq)
490 {
491  FTYPE outval;
492  FTYPE gamma,W;
493  FTYPE compute_idwmrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
494  FTYPE compute_idrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
495  FTYPE rho0,wmrho0;
496  FTYPE dwmrho0dvsq,idwmrho0dp;
497  FTYPE drho0dvsq,idrho0dp;
498  FTYPE gammasq;
499  FTYPE dpdvsq;
500 
501  gammasq=1.0+utsq;
502  gamma=sqrt(gammasq);
503  wmrho0=wmrho0_compute_utsq(Wp, D, utsq, gamma, gammasq);
504 
505  rho0=D/gamma;
506 
507  W = Wp+D;
508 
509  idwmrho0dp=compute_idwmrho0dp(whicheos,EOSextra, rho0, wmrho0);
510  dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp);
511 
512  drho0dvsq = -D*gamma*0.5; // because \rho=D/\gamma
513  idrho0dp = compute_idrho0dp(whicheos,EOSextra, rho0, wmrho0);
514 
515  dpdvsq = drho0dvsq *idrho0dp + dwmrho0dvsq *idwmrho0dp;
516 
517  // if( outval > 0. ) {
518  // dualfprintf(fail_file,"outval = %21.15g , D = %21.15g , vsq = %21.15g, W = %21.15g \n",
519  // outval, D, vsq, W );
520  // }
521 
522  return(dpdvsq);
523 }
524 
525 
527 static FTYPE dSsdvsq_calc_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE vsq)
528 {
529  FTYPE outval;
530  FTYPE gamma,W;
531  FTYPE compute_dspecificSdrho_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
532  FTYPE compute_dspecificSdwmrho0_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
533  FTYPE rho0,wmrho0;
534  FTYPE dwmrho0dvsq,idwmrho0dSs;
535  FTYPE drho0dvsq,idrho0dSs;
536  FTYPE gammasq;
537  FTYPE dSsdvsq;
538 
539 
540  gamma=1.0/sqrt(1.0-vsq);
541  gammasq=gamma*gamma;
542  wmrho0=wmrho0_compute_vsq(Wp, D, vsq, gamma, gammasq);
543 
544  rho0=D/gamma;
545 
546  // return( (gam - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / gam ) ;
547  // outval = (gam - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / gam ;
548 
549  W = Wp+D;
550 
551  // outval = (D*gamma*0.5-W)*idwmrho0dp; // old code
552 
553  idwmrho0dSs=compute_dspecificSdwmrho0_wmrho0(whicheos,EOSextra, rho0, wmrho0);
554  dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp);
555 
556  drho0dvsq = -D*gamma*0.5; // because \rho=D/\gamma
557  idrho0dSs = compute_dspecificSdrho_wmrho0(whicheos,EOSextra, rho0, wmrho0);
558 
559  dSsdvsq = drho0dvsq *idrho0dSs + dwmrho0dvsq *idwmrho0dSs;
560 
561  // if( outval > 0. ) {
562  // dualfprintf(fail_file,"outval = %21.15g , D = %21.15g , vsq = %21.15g, W = %21.15g \n",
563  // outval, D, vsq, W );
564  // }
565 
566  return(dSsdvsq);
567 }
568 
569 
571 static FTYPE dSsdvsq_calc2_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq)
572 {
573  FTYPE outval;
574  FTYPE gamma,W;
575  FTYPE compute_dspecificSdrho_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
576  FTYPE compute_dspecificSdwmrho0_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
577  FTYPE rho0,wmrho0;
578  FTYPE dwmrho0dvsq,idwmrho0dSs;
579  FTYPE drho0dvsq,idrho0dSs;
580  FTYPE gammasq;
581  FTYPE dSsdvsq;
582 
583  gammasq=1.0+utsq;
584  gamma=sqrt(gammasq);
585  wmrho0=wmrho0_compute_utsq(Wp, D, utsq, gamma, gammasq);
586 
587  rho0=D/gamma;
588 
589  W = Wp+D;
590 
591  idwmrho0dSs=compute_dspecificSdwmrho0_wmrho0(whicheos,EOSextra, rho0, wmrho0);
592  dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp);
593 
594  drho0dvsq = -D*gamma*0.5; // because \rho=D/\gamma
595  idrho0dSs = compute_dspecificSdrho_wmrho0(whicheos,EOSextra, rho0, wmrho0);
596 
597  dSsdvsq = drho0dvsq *idrho0dSs + dwmrho0dvsq *idwmrho0dSs;
598 
599  // if( outval > 0. ) {
600  // dualfprintf(fail_file,"outval = %21.15g , D = %21.15g , vsq = %21.15g, W = %21.15g \n",
601  // outval, D, vsq, W );
602  // }
603 
604  return(dSsdvsq);
605 }
606 
607 
608 #endif
609 
610 
611 
612 
613 
614 static FTYPE dpdvsq_calc_scn(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq)
615 {
616  FTYPE outval;
617  FTYPE compute_idwmrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
618  FTYPE idwmrho0dp;
619  FTYPE rho0,wmrho0;
620  FTYPE gamma,gammasq;
621  FTYPE Wp;
622 
623  Wp=W-D;
624  gamma=1.0/sqrt(1.0-vsq);
625  gammasq=gamma*gamma;
626  wmrho0=wmrho0_compute_vsq(Wp, D, vsq, gamma, gammasq);
627 
628  rho0=D/gamma;
629  idwmrho0dp=compute_idwmrho0dp(whicheos,EOSextra, rho0, wmrho0);
630 
631  // return( (gam - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / gam ) ;
632  outval = ( 0.5 * D / sqrt(1.-vsq) - W )*idwmrho0dp ;
633 
634  // if( outval > 0. ) {
635  // dualfprintf(fail_file,"outval = %21.15g , D = %21.15g , vsq = %21.15g, W = %21.15g \n",
636  // outval, D, vsq, W );
637  // }
638 
639  return(outval);
640 }
641 
644 static FTYPE wmrho0_compute_vsq(FTYPE Wp, FTYPE D, FTYPE vsq, FTYPE gamma,FTYPE gammasq)
645 {
646  return(Wp/gammasq - D*vsq/(1.0+gamma));
647 
648 }
649 
650 
653 static FTYPE wmrho0_compute_utsq(FTYPE Wp, FTYPE D, FTYPE utsq, FTYPE gamma,FTYPE gammasq)
654 {
655  return( (Wp - D*utsq/(1.0+gamma) )/gammasq );
656 
657 }
658 
659 //#undef CRAZYDEBUG
Generated on Fri May 20 2016 15:52:35 for HARM by   doxygen 1.8.3.1