Ticket #41: GrnVryDpth.c
| File GrnVryDpth.c, 1.8 kB (added by peter, 9 months ago) |
|---|
| Line | |
|---|---|
| 1 | /*GrnVryDpth set grains abundance as a function of depth into cloud*/ |
| 2 | STATIC double GrnVryDpth( |
| 3 | |
| 4 | /* nd is the number of the grain bin. The values are listed in the Cloudy output, |
| 5 | * under "Average Grain Properties", and can easily be obtained by doing a trial |
| 6 | * run without varying the grain abundance and setting stop zone to 1 */ |
| 7 | |
| 8 | long int nd) |
| 9 | { |
| 10 | double GrnVryDpth_v; |
| 11 | |
| 12 | # ifdef DEBUG_FUN |
| 13 | fputs( "<+>GrnVryDpth()\n", debug_fp ); |
| 14 | # endif |
| 15 | |
| 16 | /* set grains abundance as a function of depth into cloud |
| 17 | * NB most quantities are undefined for first calls to this sub */ |
| 18 | /* nd is the pointer to the grain bin. This routine must return |
| 19 | * a scale factor for the abundance at this position. */ |
| 20 | |
| 21 | if( gv.bin[nd]->lgDustVary ) |
| 22 | { |
| 23 | #if 0 |
| 24 | /* sample code - the scale factor will be the hydrogen ionization fraction |
| 25 | * following can be used for any other element - first dimension |
| 26 | * is atomic number (ipHYDROGEN, ipHELIUM, etc.), and second is stage |
| 27 | * of ionization (0=neutral, 1=singly ionized, etc.) - this function is |
| 28 | * probably appropriate for PAHs, which do not exist in ionized gas but |
| 29 | * are present beyond the ionization front */ |
| 30 | if( dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHYDROGEN][0] > 0. ) |
| 31 | GrnVryDpth_v = dense.xIonDense[ipHYDROGEN][0]/ |
| 32 | (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHYDROGEN][0]); |
| 33 | else |
| 34 | GrnVryDpth_v = 1.; |
| 35 | #else |
| 36 | if( phycon.te >= 4000. ) |
| 37 | GrnVryDpth_v = 1.e-10; |
| 38 | else if( phycon.te >= 3000. ) |
| 39 | GrnVryDpth_v = 1.e-10 - 1./PI*atan2((phycon.te-3000.)*(phycon.te-4000.), |
| 40 | (phycon.te-3000.)*2000.+(phycon.te-4000.)*2000.); |
| 41 | else |
| 42 | GrnVryDpth_v = 1.; |
| 43 | |
| 44 | ASSERT( GrnVryDpth_v > 0. && GrnVryDpth_v <= 1.000001 ); |
| 45 | #endif |
| 46 | } |
| 47 | else |
| 48 | { |
| 49 | GrnVryDpth_v = 1.; |
| 50 | } |
| 51 | |
| 52 | |
| 53 | # ifdef DEBUG_FUN |
| 54 | fputs( " <->GrnVryDpth()\n", debug_fp ); |
| 55 | # endif |
| 56 | return GrnVryDpth_v; |
| 57 | } |
| 58 | /*lint +e662 creation of out of bounds pointer */ |
