00001
00002
00003 #include "cddefines.h"
00004 #include "physconst.h"
00005 #include "lines_service.h"
00006 #include "elementnames.h"
00007 #include "taulines.h"
00008 #include "path.h"
00009 #include "trace.h"
00010 #include "phycon.h"
00011 #include "thermal.h"
00012 #include "dense.h"
00013 #include "iso.h"
00014 #include "helike.h"
00015 #include "helike_recom.h"
00016 #include "helike_cs.h"
00017 #include "helike_einsta.h"
00018 #include "hydroeinsta.h"
00019
00020 static double EthRyd;
00021
00022 static FILE *ioOFP;
00023
00024
00025 static char **chLevel;
00026
00027
00028
00029 static double EionWN[LIMELM] =
00030
00031
00032 {-DBL_MAX,
00033 198310.6679 ,610003.839889137,1241136.72201499,2091948.45665631,3162116.52584231,
00034 4452446.95015668,5962133.81875305,7692790.05069734,9645221.44709864,11814589.7994457,
00035 14209766.0528639,16822685.5022862,19661412.9625169,22717883.6187518,26000162.0663204,
00036 29508248.5246975,33234078.1790787,37185715.7345311,41363161.0813172,45766414.4389118,
00037 50395475.4781030,55258409.0136949,60339085.8550283,65653635.1927626,71202056.8074231,
00038 76976286.4328920,82984388.3352872,89194104.5722390,95726403.3055320};
00039
00040
00041
00042 static double EionRYD[LIMELM] =
00043
00044
00045 {-DBL_MAX,
00046 1.807387521,5.558764,11.310070,19.063237,28.815326,40.573682,54.330961,70.101861,
00047 87.893725,107.662464,129.488916,153.299590,179.167978,207.020588,236.930910,
00048 268.898946,302.851204,338.861175,376.928858,417.054255,459.237363,503.551674,
00049 549.850208,598.279945,648.840883,701.459535,756.209388,812.796486,872.323172};
00050
00051
00052
00053 #define NHE1LEVELS 111
00054
00055 static double He1Energies[NHE1LEVELS] =
00056 {0.0 , 159855.9734, 166277.4390, 169087.8298, 169086.8417, 169086.7652, 171134.8957,
00057 183236.7908, 184864.8281, 185564.6657, 186101.5615, 186104.9656, 186209.3638, 190298.6619,
00058 190940.6075, 191217.0826, 191444.4868, 191446.4547, 191451.8805, 191451.8964, 191492.7108,
00059 193346.9900, 193663.5106, 193800.7280, 193917.1538, 193918.2888, 193921.1207, 193921.1298,
00060 193921.6166, 193921.6209, 193942.4612, 194936.1184, 195114.8672, 195192.7542, 195260.0724,
00061 195260.7694, 195262.4251, 195262.4307, 195262.7236, 195262.7261, 195262.7930, 195262.7947,
00062 195274.9074, 195868.2357, 195978.8938, 196027.3216, 196069.6730, 196070.1273, 196071.1763,
00063 196071.1800, 196071.3686, 196071.3702, 196071.4141, 196071.4151, 196071.4283, 196071.4290,
00064 196079.0865, 196461.3605, 196534.5628, 196566.7159, 196595.0620, 196595.3730, 196596.0785,
00065 196596.0810, 196596.2092, 196596.2103, 196596.2404, 196596.2411, 196596.2503, 196596.2508,
00066 196596.2541, 196596.2544, 196601.3992, 196861.9861, 196912.9014, 196935.3339, 196955.2261,
00067 196955.4477, 196955.9445, 196955.9463, 196956.0373, 196956.0380, 196956.0595, 196956.0600,
00068 196956.0666, 196956.0670, 196956.0693, 196956.0696, 196956.0705, 196956.0707, 196959.6917,
00069 197145.2320, 197182.0643, 197198.3343, 197212.8252, 197212.9885, 197213.3513, 197213.3527,
00070 197213.4194, 197213.4200, 197213.4358, 197213.4362, 197213.4411, 197213.4414, 197213.4431,
00071 197213.4433, 197213.4440, 197213.4442, 197213.4445, 197213.4446, 197216.0885};
00072
00073
00074
00075 #define NIONLEVELS 31
00076
00077
00078
00079
00080
00081 static double IonEnergies[LIMELM-2][NIONLEVELS] =
00082 {
00083
00084 {0.00, 476034.98, 491374.60, 494266.57, 494261.17, 494263.44, 501808.59,
00085 554754.45, 558777.88, 559501.16, 561243.67, 561273.62, 561752.82,
00086 579981.33, 581596.77, 581886.34, 582613.64, 582630.95, 582642.97,
00087 582644.04, 582830.11, 591184.26, 591989.55, 592134.36, 592504.32,
00088 592514.43, 592520.11, 592521.11, -1.00 , -1.00 , 592634.91},
00089
00090 {0.00, 956502.00, 981178.00, 983366.00, 983355.00, 983370.00, 997454.00,
00091 1121184.00, 1127705.00, 1128300.00, 1131383.00, 1131462.00, 1132390.00,
00092 1175295.00, 1178005.00, 1178174.00, 1179451.00, 1179495.00, 1179515.00,
00093 1179514.00, 1179830.00, 1199650.00, -1.00 , 1201060.00, 1201702.00,
00094 1201800.00, 1201730.00, 1201742.00, -1.00 , -1.00 , 1201894.00},
00095
00096 {0.00, 1601545.00, 1635720.00, 1636938.00, 1636922.00, 1636975.00, 1657980.00,
00097 1882740.00, 1891790.00, 1892221.00, 1896710.00, 1896836.00, 1898063.00,
00098 1976420.00, -1.00 , 1980291.11, 1982132.67, 1982220.00, 1982262.67,
00099 1982240.00, 1982762.00, -1.00 , -1.00 , 2020730.00, 2021700.00,
00100 -1.00 , 2021665.71, 2021760.00, 2021770.00, 2021770.00, 2022044.00},
00101
00102 {0.00, 2411262.00, 2455024.00, 2455162.74, 2455150.23, 2455286.01, 2483371.00,
00103 2839562.00, 2851180.00, 2851418.00, 2857309.67, 2857529.00, 2859375.00,
00104 2983541.00, -1.00 , 2988359.00, 2990776.00, 2990923.00, 2990923.40,
00105 2990923.40, 2991710.00, 3048927.00, -1.00 , 3051332.00, 3052589.00,
00106 3052656.00, 3052653.30, 3052653.30, 3052659.40, 3052659.40, 3053044.00},
00107
00108 {0.00, 3385890.00, 3439274.00, 3438312.46, 3438321.13, 3438612.15, 3473790.00,
00109 3991860.00, -1.00 , 4006160.00, 4013460.00, 4013770.00, 4016390.00,
00110 4196800.00, 4202520.00, 4202620.00, 4205820.00, 4205830.00, 4205810.00,
00111 4205820.00, 4206810.00, 4290150.00, 4293020.00, 4293080.00, 4294570.00,
00112 4294670.00, 4294700.00, 4294700.00, -1.00 , -1.00 , 4296090.00},
00113
00114 {0.00, 4524640.00, 4588380.00, 4585620.76, 4585679.58, 4586231.19, 4629201.00,
00115 5338820.00, 5356420.00, 5355670.00, 5364422.67, 5365470.00, 5368550.00,
00116 5616340.00, 5623100.00, 5622600.00, 5626225.33, 5626670.00, 5626210.00,
00117 5626840.00, 5628100.00, 5742610.00, -1.00 , 5745440.00, 5747509.33,
00118 5748230.00, 5747200.00, 5747820.00, -1.00 , -1.00 , 5748450.00},
00119
00120 {0.00, 5830040.00, 5903100.00, 5900600.00, 5900750.00, 5901700.00, 5949900.00,
00121 6885090.00, 6903270.00, 6902560.00, 6914073.33, 6915900.00, 6916590.00,
00122 7244270.00, -1.00 , 7250390.00, 7255960.00, 7254240.00, 7256750.00,
00123 7257260.00, 7256370.00, -1.00 , -1.00 , 7410270.00, 7413940.00,
00124 7412290.00, 7414760.00, 7415300.00, 7414780.00, -1.00 , 7414000.00},
00125
00126 {0.00, 7299940.00, 7382680.00, 7378205.53, 7378506.43, 7380050.00, 7436560.00,
00127 8623000.00, 8644880.00, 8644744.44, 8657128.67, 8662400.00, 8660530.00,
00128 9075200.00, 9084060.00, 9084141.11, 9090355.33, -1.00 , 9089800.00,
00129 9094400.00, 9090630.00, 9282200.00, 9286650.00, 9286713.33, 9288500.00,
00130 -1.00 , 9289800.00, 9294400.00, -1.00 , -1.00 , 9290000.00},
00131
00132 {0.00, 8935337.00, 9027981.00, 9022354.10, 9022876.10, 9025284.70, 9088700.00,
00133 10558946.00, 10583431.00, 10583323.56, 10596783.40, 10597475.00, 10601080.00,
00134 11115065.00, 11124986.00, 11125102.78, 11130639.00, 11131017.00, 11131051.00,
00135 11131056.00, 11132393.00, 11369887.00, 11374868.00, 11374959.89, 11377767.00,
00136 11377984.00, 11377987.00, 11377991.00, -1.00 , -1.00 , 11378646.00},
00137
00138 {0.00, 10736136.00, 10838778.00, 10831985.83, 10832819.18, 10836391.13, 10906612.00,
00139 12691170.00, 12718304.00, 12718286.89, 12733392.33, 12734298.00, 12738006.00,
00140 13361991.00, 13372977.00, 13373168.22, 13379472.60, 13379830.00, 13379893.00,
00141 13379898.00, 13381265.00, 13669618.00, 13675137.00, 13675269.22, 13678467.13,
00142 13678680.00, -1.00 , -1.00 , -1.00 , -1.00 , 13679363.00},
00143
00144 {0.00, 12703061.00, 12815760.00, 12807847.00, 12809088.00, 12814213.00, 12891081.00,
00145 15020463.00, 15050257.00, 15050434.00, 15067287.07, 15068371.00, 15072141.00,
00146 15816791.00, 15828851.00, 15829158.67, 15836125.13, 15836581.00, -1.00 ,
00147 -1.00 , 15838068.00, 16182216.00, 16188281.00, 16188471.33, 16192010.13,
00148 16192244.00, -1.00 , -1.00 , -1.00 , -1.00 , 16192975.00},
00149
00150 {0.00, 14835945.00, 14958753.00, 14949756.42, 14951532.63, 14958690.57, 15042040.00,
00151 17546734.00, 17579166.00, 17579686.44, 17598406.93, 17599605.00, 17603422.00,
00152 18479389.00, 18492532.00, 18493007.56, 18500821.00, 18501245.00, -1.00 ,
00153 -1.00 , 18502736.00, 18907613.00, 18914246.00, 18914502.78, 18918476.30,
00154 18918694.00, -1.00 , -1.00 , -1.00 , -1.00 , 18919421.00},
00155
00156 {0.00, 17135768.00, 17268828.00, 17258746.00, 17261164.00, 17270908.00, 17360546.00,
00157 20271100.00, 20306284.00, 20307209.11, 20327865.87, 20329412.00, 20332952.00,
00158 21350958.00, 21365192.00, 21365892.89, 21374428.00, 21375044.00, 21375302.00,
00159 21375302.00, 21376454.00, 21846994.00, 21854144.00, 21854552.89, 21858894.80,
00160 21859210.00, 21859340.00, 21859340.00, -1.00 , -1.00 , 21859464.00},
00161
00162 {0.00, 19602076.00, 19745473.00, 19734297.61, 19737518.84, 19750576.04, 19846285.00,
00163 23193163.00, 23231087.00, 23232596.56, 23255347.53, 23257195.00, 23260416.00,
00164 24431101.00, 24446439.00, 24447429.78, 24456981.30, 24457576.00, -1.00 ,
00165 -1.00 , 24458842.00, 24999972.00, 25007605.00, 25008238.11, 25013102.40,
00166 25013407.00, -1.00 , -1.00 , -1.00 , -1.00 , 25014007.00},
00167
00168 {0.00, 22236180.00, 22390000.00, 22377820.00, 22381940.00, 22399100.00, 22500680.00,
00169 26314360.00, 26355050.00, 26357324.44, 26382328.67, 26384530.00, 26387270.00,
00170 27720900.00, 27738000.00, 27738966.67, 27749520.00, 27750400.00, 27757331.14,
00171 27757178.00, 27751600.00, 28367700.00, 28376500.00, 28376977.78, 28382286.67,
00172 28382800.00, -1.00 , -1.00 , -1.00 , -1.00 , 28383400.00},
00173
00174 {0.00, 25038230.00, 25202480.00, 25189388.10, 25194588.99, 25216810.57, 25323950.00,
00175 29634850.00, 29678210.00, 29681541.11, 29713920.00, 29715070.00, 29713740.00,
00176 31221700.00, 31239280.00, 31241087.78, -1.00 , -1.00 , -1.00 ,
00177 -1.00 , 31254280.00, 31951370.00, 31960150.00, 31961182.22, -1.00 ,
00178 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 31967860.00},
00179
00180 {0.00, 28004980.00, 28180480.00, 28165880.00, 28172670.00, 28200800.00, 28312910.00,
00181 33151930.00, 33198090.00, 33202622.22, -1.00 , -1.00 , 33237140.00,
00182 -1.00 , 34948420.00, -1.00 , -1.00 , -1.00 , -1.00 ,
00183 -1.00 , 34964920.00, 35747360.00, 35756712.00, 35758032.00, -1.00 ,
00184 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 35765192.00},
00185
00186 {0.00, 31142150.00, 31328450.00, 31312818.51, 31320486.08, 31356326.72, 31473810.00,
00187 36870940.00, 36919930.00, 36925900.00, -1.00 , -1.00 , 36962850.00,
00188 38850670.00, 38870530.00, 38873536.67, -1.00 , -1.00 , -1.00 ,
00189 -1.00 , 38888680.00, 39761380.00, 39771310.00, 39772968.89, -1.00 ,
00190 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 39780630.00},
00191
00192 {0.00, 34448120.00, 34645360.00, 34628770.00, 34638550.00, 34682810.00, 34805000.00,
00193 40790620.00, 40842480.00, 40850158.89, -1.00 , -1.00 , 40889690.00,
00194 42983370.00, 43004390.00, 43008165.56, -1.00 , -1.00 , -1.00 ,
00195 -1.00 , 43024380.00, 43992240.00, 44002740.00, 44004803.33, -1.00 ,
00196 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 44013010.00},
00197
00198 {0.00, 37923880.00, 38131640.00, 38114760.00, 38125260.00, 38180620.00, 38308340.00,
00199 44911910.00, 44966970.00, 44976472.22, 45016028.67, 45021140.00, 45018670.00,
00200 47328500.00, 47351600.00, 47355644.44, 47372326.67, 47374500.00, -1.00 ,
00201 -1.00 , 47373500.00, 48440800.00, 48452600.00, 48454581.11, -1.00 ,
00202 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 48463700.00},
00203
00204 {0.00, 41568880.00, 41787830.00, 41770130.00, 41782100.00, 41849950.00, 41982380.00,
00205 49234710.00, 49292760.00, 49304574.44, 49347922.67, 49353910.00, 49349740.00,
00206 51886600.00, 51910900.00, 51915955.56, 51934113.33, 51936800.00, -1.00 ,
00207 -1.00 , 51935100.00, 53107300.00, 53119700.00, 53122217.78, -1.00 ,
00208 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 53132000.00},
00209
00210 {0.00, 45384110.00, 45614410.00, 45595910.00, 45609360.00, 45691820.00, 45828830.00,
00211 53760100.00, 53821190.00, 53835660.00, 53883181.33, 53890160.00, 53884060.00,
00212 56658500.00, 56684100.00, 56690288.89, 56710306.67, 56713200.00, -1.00 ,
00213 -1.00 , 56710700.00, 57992700.00, 58005800.00, 58008904.44, -1.00 ,
00214 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 58019400.00},
00215
00216 {0.00, 49370240.00, 49612040.00, 49592800.00, 49607700.00, 49707130.00, 49848620.00,
00217 58488800.00, 58553000.00, 58570522.22, 58622620.00, 58630700.00, 58622500.00,
00218 61644700.00, 61671800.00, 61679233.33, 61701013.33, 61704700.00, -1.00 ,
00219 -1.00 , 61701200.00, 63097900.00, 63111600.00, 63115327.78, -1.00 ,
00220 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 63126600.00},
00221
00222 {0.00, 53527760.00, 53781230.00, 53761280.00, 53777570.00, 53896550.00, 54042490.00,
00223 63421700.00, 63489000.00, 63509966.67, 63567120.00, 63576500.00, 63565800.00,
00224 66846900.00, 66875000.00, 66884011.11, 66908120.00, 66912100.00, -1.00 ,
00225 -1.00 , 66907600.00, 68423800.00, 68438100.00, 68442730.00, -1.00 ,
00226 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 68454800.00},
00227
00228 {0.00, 57857380.00, 58122700.00, 58102090.00, 58119680.00, 58261180.00, 58411430.00,
00229 68560000.00, 68630600.00, 68655444.44, 68718160.00, 68728900.00, 68715500.00,
00230 72266000.00, 72295500.00, 72306144.44, 72332586.67, 72337100.00, -1.00 ,
00231 -1.00 , 72331500.00, 73972000.00, 73987000.00, 73992318.89, -1.00 ,
00232 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 74005400.00},
00233
00234 {0.00, 62358960.00, 62637200.00, 62615022.80, 62633778.00, 62800884.44, 62952670.00,
00235 73903340.00, 73976370.00, 74005924.44, -1.00 , -1.00 , 74070580.00,
00236 77900890.00, 77930480.00, 77943808.89, -1.00 , -1.00 , -1.00 ,
00237 -1.00 , 77970500.00, 79740940.00, 79755710.00, 79762730.00, -1.00 ,
00238 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 79776290.00},
00239
00240 {0.00, 67035380.00, 67324970.00, 67303150.00, 67322950.00, 67519170.00, 67678080.00,
00241 79453200.00, 79530300.00, 79564500.00, 79640046.67, 79654200.00, 79634300.00,
00242 83754400.00, 83786500.00, 83801155.56, 83832993.33, 83838900.00, -1.00 ,
00243 -1.00 , 83830600.00, 85734300.00, 85750700.00, 85759175.00, -1.00 ,
00244 -1.00 , -1.00 , -1.00 , -1.00 , -1.00 , 85773200.00},
00245
00246 {0.00, 71886300.00, 72188400.00, 72166200.00, 72186600.00, 72415600.00, 72579000.00,
00247 85212700.00, 85293200.00, 85332955.56, 85415826.67, 85431900.00, 85408300.00,
00248 89828600.00, 89862000.00, 89879022.22, 89913666.67, 89920800.00, -1.00 ,
00249 -1.00 , 89910900.00, 91953400.00, 91970500.00, 91979222.22, 91995073.33,
00250 92000500.00, -1.00 , -1.00 , -1.00 , -1.00 , 91994500.00}};
00251
00252
00253 static double he_energy(double Eff_n,long int nelem, long int ipLo )
00254 {
00255 double Ef;
00256
00257
00258 if( nelem == ipHELIUM && ipLo < NHE1LEVELS )
00259 {
00260 Ef = EionWN[ipHELIUM] - He1Energies[ipLo];
00261 }
00262
00263 else if( nelem > ipHELIUM && nelem <= ipZINC &&
00264 ipLo < NIONLEVELS && IonEnergies[nelem-2][ipLo] > 0. )
00265 {
00266 Ef = EionWN[nelem] - IonEnergies[nelem-2][ipLo];
00267 }
00268 else
00269
00270 {
00271
00272
00273
00274
00275 Ef = 0.999862926 * RYD_INF * (nelem/Eff_n) * (nelem/Eff_n);
00276 }
00277
00278 ASSERT(Ef > 0.);
00279 return Ef;
00280 }
00281
00282
00283 static double defect(long int nelem, long int ipLo)
00284 {
00285
00286 double qd,a,b,c;
00287
00288
00289
00290
00291
00292 double HeDefectAsymptotes[2][10] = {
00293 {1.40005E-01,-1.20673E-02,2.08056E-03,4.21484E-04,1.14868E-04,
00294 4.08648E-05,1.73548E-05,8.33891E-06,4.39680E-06,2.42075E-06},
00295 {2.97063E-01,6.81567E-02,2.82381E-03,4.27703E-04,1.17319E-04,
00296 4.25254E-05,1.85549E-05,9.24641E-06,5.30882E-06,3.02877E-06}
00297 };
00298
00299
00300
00301
00302
00303
00304 double param[3][4][3]=
00305 {
00306 {{0.6451941,0.3119437,-1.2722842},
00307 {0.7664874,0.3455675,-1.3976462},
00308 {0.8247101,0.3603131,-1.4520500},
00309 {0.8878402,0.3714450,-1.4995732}},
00310
00311 {{1.4203514,0.5311096,-2.6728087},
00312 {1.5733513,0.5997339,-2.9253834},
00313 {1.4531025,0.5924751,-2.8662756},
00314 {1.6038999,0.6342552,-3.0298071}},
00315
00316 {{-2.2323488,0.0890840,-0.5166053},
00317 {-2.0463691,0.1222081,-0.6672983},
00318 {-1.9904104,0.1328918,-0.7150879},
00319 {-1.9500974,0.1452111,-0.7649031}}
00320 };
00321
00322
00323
00324
00325
00326 double P1[4][2]=
00327 {
00328 {-56.65245,-3.661923},
00329 {-52.03411,-4.941075},
00330 {-50.43744,-5.525750},
00331 {-49.45137,-5.908615}
00332 };
00333
00334 long int n = iso.quant_desig[ipHE_LIKE][nelem][ipLo].n;
00335 long int lqn = iso.quant_desig[ipHE_LIKE][nelem][ipLo].l;
00336 long int s = iso.quant_desig[ipHE_LIKE][nelem][ipLo].s;
00337
00338 ASSERT(n >= 1L);
00339 ASSERT(n > lqn);
00340
00341 ASSERT((nelem >= ipHELIUM) && (nelem < LIMELM));
00342
00343 if( n > iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00344 {
00345
00346 qd = 0.;
00347 }
00348 else if( nelem == ipHELIUM )
00349 {
00350 if( ipLo<NHE1LEVELS && n<=iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00351 {
00352
00353 qd = n-sqrt(0.999862926*RYD_INF/(EionWN[ipHELIUM] - He1Energies[ipLo]));
00354 }
00355 else if( lqn<=9 )
00356 {
00357
00358 qd = HeDefectAsymptotes[s][lqn];
00359 }
00360 else if( s == 0 )
00361 {
00362
00363 qd = 0.0497*pow((double)lqn, -4.4303);
00364 }
00365 else
00366 {
00367
00368 qd = 0.0656*pow((double)lqn, -4.5606);
00369 }
00370 }
00371 else if( ipLo == ipHe1s1S )
00372 {
00373
00374
00375
00376 ASSERT(nelem>ipHYDROGEN && nelem<LIMELM );
00377 qd = 1.0 - nelem * sqrt(1/EionRYD[nelem]);
00378 }
00379 else
00380 {
00381
00382
00383 if( n > 5L )
00384 {
00385 n = 5L;
00386 }
00387
00388 if( lqn==1L && s==0L )
00389 {
00390 qd = 1./(P1[n-2][0] + P1[n-2][1] * (nelem+1) * log((double)nelem+1.) );
00391 }
00392
00393 else if( lqn < 2L )
00394 {
00395 a = param[2*(lqn+1)-s-1][n-2][0];
00396 b = param[2*(lqn+1)-s-1][n-2][1];
00397 c = param[2*(lqn+1)-s-1][n-2][2];
00398 qd = exp((a+c*(nelem+1))/(1.0+b*(nelem+1)));
00399 }
00400
00401
00402
00403
00404
00405
00406
00407 else
00408 {
00409 ASSERT( lqn >= 2L );
00410 qd = ( ( 0.0612/(double)nelem ) / pow((double)lqn, 4.44) );
00411 }
00412 }
00413
00414 return qd;
00415 }
00416
00417
00418
00419 static void he_assign( long int nelem)
00420 {
00421
00422
00423 long int in, il, is, ij, i = 0 , level;
00424
00425
00426 for( in = 1L; in <= iso.n_HighestResolved_max[ipHE_LIKE][nelem]; ++in )
00427 {
00428 for( il = 0L; il < in; ++il )
00429 {
00430 for( is = 1L; is >= 0; --is )
00431 {
00432
00433
00434
00435 if( (il == 1L) && (is == 0L) )
00436 continue;
00437
00438 if( (in == 1L) && (is == 1L) )
00439 continue;
00440 ij = 0;
00441 do
00442 {
00443 iso.quant_desig[ipHE_LIKE][nelem][i].n = in;
00444 iso.quant_desig[ipHE_LIKE][nelem][i].s = is;
00445 iso.quant_desig[ipHE_LIKE][nelem][i].l = il;
00446 QuantumNumbers2Index[nelem][in][il][is] = i;
00447 ++i;
00448 ++ij;
00449
00450 } while ( (in == 2) && (il == 1) && (is == 1) && (ij < 3) );
00451 }
00452 }
00453
00454 if( in > 1L )
00455 {
00456 iso.quant_desig[ipHE_LIKE][nelem][i].n = in;
00457 iso.quant_desig[ipHE_LIKE][nelem][i].s = 0L;
00458 iso.quant_desig[ipHE_LIKE][nelem][i].l = 1L;
00459 QuantumNumbers2Index[nelem][in][1][0] = i;
00460 ++i;
00461 }
00462 }
00463
00464 in = iso.n_HighestResolved_max[ipHE_LIKE][nelem] + 1;
00465 for( level = i; level< iso.numLevels_max[ipHE_LIKE][nelem]; ++level)
00466 {
00467 iso.quant_desig[ipHE_LIKE][nelem][level].n = in;
00468 iso.quant_desig[ipHE_LIKE][nelem][level].s = -LONG_MAX;
00469 iso.quant_desig[ipHE_LIKE][nelem][level].l = -LONG_MAX;
00470
00471 for( il = 0; il < in; ++il )
00472 {
00473 for( is = 0; is <= 1; ++is )
00474 {
00475 QuantumNumbers2Index[nelem][in][il][is] = level;
00476 }
00477 }
00478 ++in;
00479 }
00480 --in;
00481
00482
00483 ASSERT( i <= iso.numLevels_max[ipHE_LIKE][nelem] );
00484
00485
00486 ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) ) );
00487
00488
00489 for( in = 2L; in <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]; ++in )
00490 {
00491 for( il = 0L; il < in; ++il )
00492 {
00493 for( is = 1L; is >= 0; --is )
00494 {
00495
00496 if( (in == 1L) && (is == 1L) )
00497 continue;
00498
00499 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].n == in );
00500 if( in <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00501 {
00502
00503
00504 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].l == il );
00505 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].s == is );
00506 }
00507 }
00508 }
00509 }
00510
00511 return;
00512 }
00513
00514
00515
00516 static void printCustomAs( void )
00517 {
00518 long int i, ipLo, familyIndex, nHi, upperIndex, nelem;
00519 nelem = ipHELIUM;
00520
00521 if( (ioOFP = fopen("CustomAs.txt","w")) == NULL )
00522 return;
00523 fprintf(ioOFP,"ipLo\tnLo\tlLo\tsLo\tlHi\t" );
00524 for( i = 2; i<= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM]; i++ )
00525 {
00526 fprintf( ioOFP,"%li\t", i);
00527 }
00528 fprintf( ioOFP, "\n" );
00529
00530 for( ipLo = 0; ipLo< iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] - 1; ipLo++ )
00531 {
00532 if( N_(ipLo) < iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] )
00533 {
00534 for( familyIndex = N_(ipLo)*(N_(ipLo)+1)+1; familyIndex < (N_(ipLo)+1)*(N_(ipLo)+2)+1; familyIndex++ )
00535 {
00536 if( N_(familyIndex) <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] )
00537 {
00538 if( ( abs(L_(ipLo)-L_(familyIndex)) == 1 ) && (S_(ipLo) == S_(familyIndex)) )
00539 {
00540 fprintf( ioOFP, "%li\t%li\t%li\t%li\t%li\t",
00541 ipLo, N_(ipLo), L_(ipLo), S_(ipLo), L_(familyIndex) );
00542
00543 for( nHi = 2; nHi < N_(familyIndex); nHi++ )
00544 {
00545 fprintf( ioOFP, "\t" );
00546 }
00547
00548 for( nHi = N_(familyIndex); nHi <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM]; nHi++ )
00549 {
00550 upperIndex = QuantumNumbers2Index[ipHELIUM][nHi][L_(familyIndex)][S_(familyIndex)];
00551
00552 fprintf( ioOFP, "%.4e\t", EmisLines[ipHE_LIKE][nelem][upperIndex][ipLo].Aul );
00553 }
00554 fprintf( ioOFP, "\n" );
00555 }
00556 }
00557 }
00558 }
00559 }
00560 }
00561
00562
00563 void HeCreate(void)
00564 {
00565 double **energies, *n_effs, **SumAPerN, **RobbinsC, ***HeliumAs;
00566
00567 long int i, i1;
00568 long int j, ipLo, ipHi, ipFirstCollapsed, nelem;
00569
00570 static int nCalled = 0;
00571
00572 # define chLine_LENGTH 1000
00573 char chLine[chLine_LENGTH] ,
00574
00575 chFilename[FILENAME_PATH_LENGTH_2];
00576
00577 FILE *ioDATA;
00578 bool lgEOL, lgHIT;
00579
00580 char chSpin[2]={'1','3'};
00581 char chL[6]={'S','P','D','F','G','H'};
00582
00583 DEBUG_ENTRY( "HeCreate()" );
00584
00585
00586
00587 if( nCalled > 0 )
00588 {
00589 DEBUG_EXIT( "HeCreate()" );
00590 return;
00591 }
00592 ++nCalled;
00593
00594
00595 if( helike.lgCompileRecomb )
00596 {
00597 helike.lgNoRecombInterp = false;
00598 }
00599
00600
00601 if( helike.lgSetBenjamin )
00602 {
00603 max_num_levels = iso.numLevels_max[ipHE_LIKE][ipHELIUM];
00604 max_n = iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM];
00605 iso.lgInd2nu_On = false;
00606
00607 if( helike.lgCompileRecomb )
00608 {
00609 fprintf( ioQQQ, " Not compiling He recomb file...can not do with Benjamin command set.");
00610 helike.lgCompileRecomb = false;
00611 }
00612 if( helike.lgFSM )
00613 {
00614 fprintf( ioQQQ, " Can not include f-s mixing with the Benjamin command! I have turned it off.");
00615 helike.lgFSM = false;
00616 }
00617 }
00618 else
00619 {
00620
00621 max_num_levels = 0;
00622 max_n = 0;
00623 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00624 {
00625
00626 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00627 {
00628 max_num_levels = MAX2( max_num_levels, iso.numLevels_max[ipHE_LIKE][nelem] );
00629 max_n = MAX2( max_n, iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] );
00630 }
00631 }
00632 }
00633
00634
00635 if( helike.lgRandErrGen )
00636 {
00637
00638
00639
00640
00641
00642
00643
00644 helike.Error = (float****)MALLOC(sizeof(float***)*(unsigned)LIMELM );
00645
00646 helike.ErrorFactor = (float****)MALLOC(sizeof(float***)*(unsigned)LIMELM );
00647 }
00648
00649
00650 helike.qTot2TripS = (double*)MALLOC(sizeof(double)*(unsigned)LIMELM );
00651
00652
00653
00654 helike.cs_proton = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM );
00655 helike.cs_heplus = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM );
00656
00657
00658
00659 #if 0
00660
00661
00662
00663 helike.cs_elec_power = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM );
00664 helike.cs_prot_power = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM );
00665 helike.cs_heplus_power = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM );
00666 #endif
00667
00668
00669 HeliumAs = (double***)MALLOC(sizeof(double**)*(unsigned)(LIMELM) );
00670
00671
00672 helike.BranchRatio = (double***)MALLOC(sizeof(double**)*(unsigned)(LIMELM) );
00673
00674
00675
00676 helike.CascadeProb = (double***)MALLOC(sizeof(double**)*(unsigned)(LIMELM) );
00677
00678
00679 chLevel = (char **)MALLOC(sizeof(char *)*(unsigned)(max_num_levels) );
00680
00681
00682 RobbinsC = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) );
00683
00684
00685
00686 helike.RadEffec = (double*)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
00687
00688 if( helike.lgRandErrGen && helike.lgHugeCaseB )
00689 {
00690
00691
00692 helike.SigmaCascadeProb = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) );
00693 helike.SigmaRadEffec = (double*)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
00694 }
00695
00696
00697 energies = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
00698
00699
00700 helike.Lifetime = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
00701
00702 if( helike.lgRandErrGen && helike.lgHugeCaseB )
00703 {
00704
00705 helike.SigmaAtot = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
00706 }
00707
00708
00709 SumAPerN = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
00710
00711
00712
00713 QuantumNumbers2Index = (long****)MALLOC(sizeof(long***)*(unsigned)LIMELM );
00714
00715 for( ipLo=ipHe1s1S; ipLo < iso.numLevels_max[ipHE_LIKE][ipHELIUM];++ipLo )
00716 {
00717 chLevel[ipLo] = (char*)MALLOC(sizeof(char)*(unsigned)10 );
00718 }
00719
00720
00721
00722 for( ipLo=ipHe1s1S; ipLo < max_num_levels;++ipLo )
00723 {
00724 RobbinsC[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) );
00725
00726 if( helike.lgRandErrGen && helike.lgHugeCaseB )
00727 {
00728
00729 helike.SigmaCascadeProb[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) );
00730 }
00731 }
00732
00733 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
00734 {
00735 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00736 {
00737
00738 energies[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) );
00739
00740
00741
00742 helike.CascadeProb[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) );
00743
00744
00745 if( helike.lgRandErrGen )
00746 {
00747
00748 helike.Error[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+2) );
00749
00750 helike.ErrorFactor[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+2) );
00751
00752 for( i = 1; i< iso.numLevels_max[ipHE_LIKE][nelem] + 1; i++ )
00753 {
00754 helike.Error[nelem][i] = (float**)MALLOC(sizeof(float*)*(unsigned)(i+1) );
00755
00756 helike.ErrorFactor[nelem][i] = (float**)MALLOC(sizeof(float*)*(unsigned)(i+1) );
00757
00758 for( j = 0; j<i+1; j++ )
00759 {
00760 helike.Error[nelem][i][j] = (float*)MALLOC(sizeof(float)*(unsigned)(3) );
00761
00762 helike.ErrorFactor[nelem][i][j] = (float*)MALLOC(sizeof(float)*(unsigned)(3) );
00763
00764
00765 for( i1=0; i1<3; i1++ )
00766 {
00767 helike.Error[nelem][i][j][i1] = -1.f;
00768 helike.ErrorFactor[nelem][i][j][i1] = -1.f;
00769 }
00770 }
00771 }
00772
00773 if( helike.lgHugeCaseB )
00774 {
00775 helike.SigmaAtot[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) );
00776 }
00777 }
00778
00779
00780 helike.cs_proton[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) );
00781
00782 helike.cs_heplus[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) );
00783
00784
00785
00786 #if 0
00787 helike.cs_elec_power[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) );
00788
00789 helike.cs_prot_power[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) );
00790
00791 helike.cs_heplus_power[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) );
00792 #endif
00793
00794 HeliumAs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) );
00795
00796 helike.BranchRatio[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) );
00797
00798 for( i = 1; i< iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
00799 {
00800
00801 helike.cs_proton[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00802
00803 helike.cs_heplus[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00804
00805
00806 #if 0
00807 helike.cs_elec_power[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00808
00809 helike.cs_prot_power[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00810
00811 helike.cs_heplus_power[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00812 #endif
00813
00814 HeliumAs[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00815
00816 helike.BranchRatio[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00817
00818 for( j = 0; j< i; j++ )
00819 {
00820 helike.cs_proton[nelem][i][j] = -FLT_MAX;
00821 helike.cs_heplus[nelem][i][j] = -FLT_MAX;
00822 #if 0
00823 helike.cs_elec_power[nelem][i][j] = -FLT_MAX;
00824 helike.cs_prot_power[nelem][i][j] = -FLT_MAX;
00825 helike.cs_heplus_power[nelem][i][j] = -FLT_MAX;
00826 #endif
00827 HeliumAs[nelem][i][j] = -DBL_MAX;
00828 helike.BranchRatio[nelem][i][j] = -DBL_MAX;
00829 }
00830 }
00831
00832
00833 helike.Lifetime[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) );
00834
00835
00836 SumAPerN[nelem] = (double*)MALLOC(sizeof(double)* (unsigned)(iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) );
00837
00838
00839 QuantumNumbers2Index[nelem] = (long***)MALLOC(sizeof(long**)* (unsigned)(iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) );
00840
00841
00842 for( i = 0; i < iso.numLevels_max[ipHE_LIKE][nelem]; ++i)
00843 {
00844 helike.CascadeProb[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) );
00845
00846 for( j = 0; j< i; j++ )
00847 {
00848 helike.CascadeProb[nelem][i][j] = -DBL_MAX;
00849 }
00850 }
00851
00852 for( i = 0; i <= (iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]); ++i)
00853 {
00854
00855 QuantumNumbers2Index[nelem][i] = (long**)MALLOC(sizeof(long*)*(unsigned)MAX2(1,i) );
00856
00857 for( i1 = 0; i1 < i; ++i1 )
00858 {
00859
00860 QuantumNumbers2Index[nelem][i][i1] = (long*)MALLOC(sizeof(long)*(unsigned)2 );
00861 }
00862 }
00863 }
00864 }
00865
00866 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00867 {
00868
00869 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00870 {
00871
00872
00873 he_assign(nelem);
00874 }
00875 }
00876
00877
00878
00879 n_effs= (double *)MALLOC((unsigned)(max_num_levels+1)*sizeof(double));
00880 {
00881
00882
00883
00884 enum {DEBUG_LOC=false};
00885
00886 if( DEBUG_LOC )
00887 {
00888 for( nelem = ipHELIUM; nelem <LIMELM; ++nelem )
00889 {
00890
00891 for( i=0; i < iso.numLevels_max[ipHE_LIKE][nelem]; ++i )
00892 {
00893 fprintf(ioQQQ,"nelem %li i %li n %li s %li l %li\n",
00894 nelem, i, N_(i), S_(i), L_(i) );
00895 }
00896 if( DEBUG_LOC )
00897 cdEXIT(EXIT_SUCCESS);
00898 }
00899 }
00900 }
00901
00902
00903 HelikeTransProbSetup();
00904
00905
00906 for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
00907 {
00908
00909 double z4 = POW2((double)nelem);
00910 z4 *= z4;
00911
00912
00913 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00914 {
00915 for( i=0; i <= ( iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] ); ++i)
00916 {
00917 SumAPerN[nelem][i] = 0.;
00918 }
00919
00920
00921 ipFirstCollapsed = iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem];
00922
00923
00924 RobbinsC[ipHe1s1S][ipHe1s1S] = 1.;
00925 helike.Lifetime[nelem][ipHe1s1S] = 0.;
00926 helike.CascadeProb[nelem][ipHe1s1S][ipHe1s1S] = 1.;
00927 if( helike.lgRandErrGen && helike.lgHugeCaseB )
00928 {
00929 helike.SigmaAtot[nelem][ipHe1s1S] = 0.;
00930 helike.SigmaCascadeProb[ipHe1s1S][ipHe1s1S] = 0.;
00931 }
00932 iso.stat[ipHE_LIKE][nelem][ipHe1s1S] = 1.f;
00933 n_effs[ipHe1s1S] = N_(ipHe1s1S) - defect(nelem,ipHe1s1S);
00934 energies[nelem][ipHe1s1S] = EionWN[nelem];
00935 iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHe1s1S] = energies[nelem][ipHe1s1S]*WAVNRYD;
00936 ASSERT( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHe1s1S] > 0. );
00937
00938
00941
00942 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00943 {
00944
00945 n_effs[ipHi] = N_(ipHi) - defect(nelem,ipHi);
00946 ASSERT( ( L_(ipHi)==1 && S_(ipHi)==0 ) ||
00947 ( N_(ipHi) - n_effs[ipHi] >= 0. ) );
00948
00949 HeliumAs[nelem][ipHi][ipHi] = 0.;
00950
00951 if(ipHi >= ipFirstCollapsed )
00952 {
00953
00954 iso.stat[ipHE_LIKE][nelem][ipHi] = 4.f*POW2((float)N_(ipHi));
00955 }
00956 else if( ipHi>=ipHe2p3P0 && ipHi<=ipHe2p3P2 )
00957 {
00958
00959 iso.stat[ipHE_LIKE][nelem][ipHi] = (2.f*J_(ipHi) + 1.f);
00960 }
00961 else
00962 {
00963
00964 iso.stat[ipHE_LIKE][nelem][ipHi] = (2.f*L_(ipHi) + 1.f) * (2.f*S_(ipHi) + 1.f);
00965 }
00966
00967 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
00968 {
00969 double Enerwn, Aul, Aul1;
00970
00971 if( helike.lgHugeCaseB == false )
00972 {
00973
00974 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].nelem = (int)(nelem+1);
00975 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].IonStg = (int)(nelem);
00976
00977
00978 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo = iso.stat[ipHE_LIKE][nelem][ipLo];
00979 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi = iso.stat[ipHE_LIKE][nelem][ipHi];
00980 ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi> 0.);
00981 ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo> 0.);
00982 }
00983
00984
00985
00986
00987 {
00988
00989 if( nelem==ipHELIUM && ipHi<NHE1LEVELS &&
00990 N_(ipHi)<=iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00991 {
00992
00993
00994 energies[nelem][ipHi] = EionWN[ipHELIUM] - He1Energies[ipHi];
00995 }
00996 else if( N_(ipHi) > iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00997 {
00998
00999
01000
01001
01002 energies[nelem][ipHi] = 0.999862926* RYD_INF * POW2((double)nelem/(double)N_(ipHi));
01003 }
01004 else
01005 {
01006
01007 energies[nelem][ipHi] = he_energy(n_effs[ipHi],nelem,ipHi);
01008 }
01009
01010
01011 EthRyd = energies[nelem][ipHi] * WAVNRYD;
01012
01013
01014 iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi] = EthRyd;
01015 ASSERT( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi] > 0. );
01016
01017
01018
01019 Enerwn = energies[nelem][ipLo] - energies[nelem][ipHi];
01020
01021
01022
01023 if( Enerwn == 0 )
01024 Enerwn = 0.0001;
01025
01026 if( Enerwn < 0. )
01027 Enerwn = -1.0 * Enerwn;
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039 if( helike.lgHugeCaseB == false )
01040 {
01041 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN = (float)Enerwn;
01042 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg = (float)(Enerwn*WAVNRYD*EN1RYD);
01043 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK = (float)(Enerwn*WAVNRYD*TE1RYD );
01044 }
01045 }
01046
01047
01048
01049
01050 {
01051 if( ipHi >= ipFirstCollapsed )
01052 {
01053 if( ipLo >= ipFirstCollapsed )
01054 {
01055
01056 Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
01057 putError(nelem,ipHi,ipLo,IPRAD,0.001f);
01058
01059 ASSERT( Aul > 0.);
01060 }
01061 else
01062 {
01063
01064
01065 Aul = he_1trans( nelem, Enerwn, n_effs[ipHi], L_(ipLo)+1,
01066 S_(ipLo), -1, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
01067 Aul *= (2.*L_(ipLo)+3.) * (2.*S_(ipLo)+1.) /
01068 (4.*(double)N_(ipHi)*(double)N_(ipHi));
01069 if( L_(ipLo) != 0 )
01070 {
01071
01072 Aul1 = he_1trans( nelem, Enerwn, n_effs[ipHi], L_(ipLo)-1,
01073 S_(ipLo), -1, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
01074
01075 Aul += Aul1*(2.*L_(ipLo)-1.) * (2.*S_(ipLo)+1.) /
01076 (4.*(double)N_(ipHi)*(double)N_(ipHi));
01077 }
01078 putError(nelem,ipHi,ipLo,IPRAD,0.01f);
01079 ASSERT( Aul > 0.);
01080 }
01081
01082 if( N_(ipHi) > N_(ipLo) )
01083 SumAPerN[nelem][N_(ipHi)] += Aul;
01084 }
01085 else
01086 {
01087
01088 if( Enerwn < 0. )
01089 {
01090 Aul = he_1trans( nelem, -1.*Enerwn, n_effs[ipLo],
01091 L_(ipLo), S_(ipLo), ipLo-3, n_effs[ipHi], L_(ipHi), S_(ipHi), ipHi-3);
01092 }
01093 else
01094 {
01095 Aul = he_1trans( nelem, Enerwn, n_effs[ipHi],
01096 L_(ipHi), S_(ipHi), ipHi-3, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
01097 }
01098
01099 if( N_(ipHi) > N_(ipLo) )
01100 {
01101
01102
01103
01104 SumAPerN[nelem][N_(ipHi)] += Aul* iso.stat[ipHE_LIKE][nelem][ipHi]/
01105 (4. * (double)N_(ipHi) * (double)N_(ipHi));
01106 }
01107 }
01108
01109
01110 if( helike.lgHugeCaseB == false )
01111 {
01112
01113 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul = (float)Aul;
01114
01115 ASSERT(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul > 0.);
01116 }
01117
01118 HeliumAs[nelem][ipHi][ipLo] = Aul;
01119 ASSERT(HeliumAs[nelem][ipHi][ipLo] > 0.);
01120 }
01121 }
01122 }
01123
01124
01125
01126
01127 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
01128 {
01129 helike.Lifetime[nelem][ipHi] = 0.;
01130 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
01131 {
01132
01133
01134
01135 if( helike.lgHugeCaseB == false )
01136 {
01137 helike.Lifetime[nelem][ipHi] += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul;
01138 }
01139 else
01140 {
01141 helike.Lifetime[nelem][ipHi] += HeliumAs[nelem][ipHi][ipLo];
01142 }
01143 }
01144
01145
01146 helike.Lifetime[nelem][ipHi] = 1./helike.Lifetime[nelem][ipHi];
01147 }
01148
01149
01150
01151
01152 if( helike.lgFSM == true && helike.lgHugeCaseB == false )
01153 {
01154 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
01155 {
01156 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
01157 {
01158 DoFSMixing( nelem, ipLo, ipHi );
01159 }
01160 }
01161 }
01162
01163
01164
01165
01166 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
01167 {
01168
01169 RobbinsC[ipHi][ipHi] = 1.;
01170 helike.CascadeProb[nelem][ipHi][ipHi] = 1.;
01171 if( helike.lgRandErrGen && helike.lgHugeCaseB )
01172 {
01173 helike.SigmaAtot[nelem][ipHi] = 0.;
01174 helike.SigmaCascadeProb[ipHi][ipHi] = 0.;
01175 }
01176
01177 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
01178 {
01179 RobbinsC[ipHi][ipLo] = 0.;
01180 helike.CascadeProb[nelem][ipHi][ipLo] = 0.;
01181
01182 if( helike.lgHugeCaseB == false )
01183 {
01184
01185
01186
01187 ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN > 0. ||
01188 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul <= iso.SmallA );
01189
01190
01191 if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul <= iso.SmallA ||
01192 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN <= 0.)
01193 {
01194
01195
01196
01197
01198
01199 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng = 1e6f;
01200
01201 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf = 1e-20f;
01202
01203
01204 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity = 0.;
01205 }
01206 else
01207 {
01208
01209
01210 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng =
01211 (float)fabs(1.0e8/
01212 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN/
01213 RefIndex( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN ));
01214
01215 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf =
01216 (float)(GetGF(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul,
01217 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
01218 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi));
01219 ASSERT(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf > 0.);
01220
01221
01222 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity =
01223 (float)(abscf(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf,
01224 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
01225 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo));
01226 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity =
01227 MAX2(0.f, EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity );
01228 }
01229 {
01230
01231
01232
01233
01234 enum {DEBUG_LOC=false};
01235
01236 if( DEBUG_LOC )
01237 {
01238 if( (nelem == ipHELIUM) && (EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul> 1e-19) &&
01239 N_(ipLo) <= 5 && N_(ipHi)<=10 )
01240 {
01241 fprintf(ioQQQ,"Z %li Lo %li n %li s %li l %li \t ",
01242 nelem+1, ipLo, N_(ipLo), S_(ipLo), L_(ipLo) );
01243 fprintf(ioQQQ," Hi %li n %li s %li l %li \t",
01244 ipHi, N_(ipHi), S_(ipHi), L_(ipHi) );
01245 fprintf(ioQQQ,"%.4e\t%.4e\tcs\t%.4e\n",
01246
01247 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
01248 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul ,
01249 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs);
01250 }
01251 }
01252 }
01253
01254
01255 if( ipLo == ipHe1s1S )
01256 {
01257 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = ipPRD;
01258 }
01259 else
01260 {
01261
01262
01263
01264 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = ipCRDW;
01265 }
01266
01267 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ipCont = INT_MIN;
01268 }
01269
01270
01271 iso.Pop2Ion[ipHE_LIKE][nelem][ipLo] = 0.;
01272
01273 if( helike.lgRandErrGen && helike.lgHugeCaseB )
01274 {
01275 ASSERT( helike.Error[nelem][ipHi][ipLo][IPRAD] >= 0. );
01276
01277 helike.SigmaAtot[nelem][ipHi] +=
01278 pow( helike.Error[nelem][ipHi][ipLo][IPRAD] * HeliumAs[nelem][ipHi][ipLo], 2. );
01279 }
01280 }
01281
01282 helike.BranchRatio[nelem][ipHi][ipHe1s1S] = 0.;
01283
01284 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
01285 {
01286 helike.BranchRatio[nelem][ipHi][ipLo] = HeliumAs[nelem][ipHi][ipLo] * helike.Lifetime[nelem][ipHi];
01287
01288
01289
01290
01291
01292
01293
01294
01295 ASSERT( helike.BranchRatio[nelem][ipHi][ipLo] <= 1.007 );
01296 }
01297
01298 if( helike.lgRandErrGen && helike.lgHugeCaseB )
01299 {
01300
01301 helike.SigmaAtot[nelem][ipHi] = sqrt( helike.SigmaAtot[nelem][ipHi] );
01302 }
01303
01304
01305
01306 {
01307 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
01308 {
01309 double SigmaA, SigmaCul = 0.;
01310
01311
01312 for( i=ipLo ; i<ipHi; i++ )
01313 {
01314 RobbinsC[ipHi][ipLo] += helike.BranchRatio[nelem][ipHi][i] * RobbinsC[i][ipLo];
01315 if( helike.lgRandErrGen && helike.lgHugeCaseB )
01316 {
01317
01318 SigmaA = helike.Error[nelem][ipHi][i][IPRAD]*HeliumAs[nelem][ipHi][i];
01319 SigmaCul +=
01320 pow(SigmaA*RobbinsC[i][ipLo]*helike.Lifetime[nelem][ipHi], 2.) +
01321 pow(helike.SigmaAtot[nelem][ipHi]*helike.BranchRatio[nelem][ipHi][i]*RobbinsC[i][ipLo]*helike.Lifetime[nelem][ipHi], 2.) +
01322 pow(helike.SigmaCascadeProb[i][ipLo]*helike.BranchRatio[nelem][ipHi][i], 2.);
01323 }
01324
01325 }
01326 helike.CascadeProb[nelem][ipHi][ipLo] = RobbinsC[ipHi][ipLo];
01327 if( helike.lgRandErrGen && helike.lgHugeCaseB )
01328 {
01329
01330 SigmaCul = sqrt(SigmaCul);
01331 helike.SigmaCascadeProb[ipHi][ipLo] = SigmaCul;
01332 }
01333
01334 }
01335 }
01336 }
01337
01338
01342
01343 {
01344
01345 enum {DEBUG_LOC=false};
01346
01347
01348 if( DEBUG_LOC && (nelem == ipHELIUM ) )
01349 {
01350
01351 long int hi_l,hi_s;
01352 double Bm;
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365 hi_l = 1;
01366 hi_s = 1;
01367 ipLo = ipHe2s3S;
01368
01369 ASSERT( hi_l != iso.quant_desig[ipHE_LIKE][nelem][ipLo].l );
01370
01371 fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
01372 fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");
01373
01374 for( ipHi=ipHe2p3P2; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; ipHi++ )
01375 {
01376
01377 if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].l == 1) && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == 1) )
01378 {
01379 fprintf(ioQQQ,"\n%ld\t",iso.quant_desig[ipHE_LIKE][nelem][ipHi].n);
01380 j = 0;
01381 Bm = 0;
01382 for( i = ipLo; i<=ipHi; i++)
01383 {
01384 if( (iso.quant_desig[ipHE_LIKE][nelem][i].l == hi_l) && (iso.quant_desig[ipHE_LIKE][nelem][i].s == hi_s) )
01385 {
01386 if( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) )
01387 {
01388 Bm += RobbinsC[ipHi][i] * ( helike.BranchRatio[nelem][i][ipHe2p3P0] +
01389 helike.BranchRatio[nelem][i][ipHe2p3P1] + helike.BranchRatio[nelem][i][ipHe2p3P2] );
01390 }
01391 else
01392 Bm += RobbinsC[ipHi][i] * helike.BranchRatio[nelem][i][ipLo];
01393
01394 if( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) )
01395 {
01396 j++;
01397 if(j == 3)
01398 {
01399 fprintf(ioQQQ,"%2.4e\t",Bm);
01400 Bm = 0;
01401 }
01402 }
01403 else
01404 {
01405 fprintf(ioQQQ,"%2.4e\t",Bm);
01406 Bm = 0;
01407 }
01408 }
01409 }
01410 }
01411 }
01412 fprintf(ioQQQ,"\n\n");
01413 }
01414 }
01415
01416
01417
01418
01419 ipLo = ipHe1s1S;
01420
01421 for( ipHi=2; ipHi <iso.nLyman[ipHE_LIKE]; ipHi++ )
01422 {
01423 double Enerwn , Aul ;
01424 float gHi , gLo;
01425
01426
01427
01428
01429 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].nelem = (int)(nelem+1);
01430 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].IonStg = (int)(nelem);
01431
01432
01433 gHi = 3.f;
01434
01435 gLo = 1.f;
01436
01437 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gLo = gLo;
01438 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gHi = gHi;
01439
01440
01441
01442 Enerwn = EionWN[nelem] - 0.999862926*RYD_INF*POW2((double)nelem)/POW2((double)ipHi);
01443
01444 if( nelem == ipHELIUM )
01445 {
01446
01447 Aul = (1.508e10) / pow((double)ipHi,2.975);
01448 }
01449 else
01450 {
01451
01452
01453
01454
01455 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)ipHi,3.1);
01456
01457 }
01458
01459 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul = (float)Aul;
01460
01461
01462
01463 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN = (float)Enerwn;
01464
01465 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyErg =
01466 (float)(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN * WAVNRYD*
01467 EN1RYD);
01468
01469
01470 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].dampXvel = (float)(Aul/
01471 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN/PI4);
01472
01473 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyK =
01474 (float)(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN * WAVNRYD*
01475 TE1RYD );
01476
01477
01478
01479
01480
01481 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].WLAng =
01482 (float)fabs(1.0e8/
01483 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN/
01484 RefIndex( iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN));
01485
01486 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gf =
01487 (float)(GetGF(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul,
01488 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN,
01489 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gHi));
01490
01491
01492 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].opacity =
01493 (float)(abscf(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gf,
01494 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN,
01495 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gLo));
01496
01497 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].iRedisFun = iso.ipResoRedist[ipHE_LIKE];
01498
01499
01500 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pesc = 1.0;
01501 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pelec_esc = 0.0;
01502
01503
01504 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pdest = 0.0;
01505
01506
01507 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].PopHi = 0.0;
01508 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].PopLo = 0.0;
01509
01510
01511 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].xIntensity = 0.0;
01512
01513
01514
01515 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].pump = 0.;
01516 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].cool = 0.;
01517 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].heat = 0.;
01518 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].ColOvTot = 0.;
01519 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].PopOpc = 0.;
01520 {
01521
01522
01523
01524
01525 enum {DEBUG_LOC=false};
01526
01527 if( DEBUG_LOC )
01528 {
01529 fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
01530 nelem+1,
01531 ipHi,
01532 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul ,
01533 iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].opacity
01534 );
01535 }
01536 }
01537 }
01538
01539
01540
01541 {
01542
01543 enum {DEBUG_LOC=false};
01544
01545 if( DEBUG_LOC && (nelem == ipHELIUM) )
01546 {
01547 printCustomAs();
01548 }
01549 }
01550
01551
01552
01553
01554
01555 {
01556 for( i=2; i < iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] - 1; ++i)
01557 {
01558 ASSERT( SumAPerN[nelem][i] > SumAPerN[nelem][i+1] );
01559 }
01560 }
01561 {
01562
01563 enum {DEBUG_LOC=false};
01564
01565 if( DEBUG_LOC )
01566 {
01567 for( i = 2; i<= (iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]); ++i)
01568 {
01569 fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[nelem][i]);
01570 }
01571 }
01572 }
01573
01574
01575 if( helike.lgHugeCaseB == false )
01576 {
01577
01578 EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].WLAng = 0.;
01579
01580
01583 EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].opacity /= 1e4f;
01584
01585
01586 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
01587 {
01588 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
01589 {
01590
01591 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].dampXvel = (float)(
01592 1./helike.Lifetime[nelem][ipHi]/
01593 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN/PI4);
01594 ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].dampXvel >= 0. );
01595
01596
01597
01598 }
01599 }
01600 }
01601 }
01602 }
01603
01604
01605
01606
01607 HelikeRecombSetup();
01608
01609
01610
01611
01612 free(n_effs );
01613
01614
01615 for( ipLo=ipHe1s1S; ipLo < iso.numLevels_max[ipHE_LIKE][ipHELIUM]-iso.nCollapsed_max[ipHE_LIKE][ipHELIUM]; ++ipLo )
01616 {
01617 nelem = ipHELIUM;
01618 sprintf( chLevel[ipLo] , "%li %c%c", N_(ipLo), chSpin[S_(ipLo)], chL[MIN2(5,L_(ipLo))] );
01619 }
01620
01621
01622 if( (trace.lgTrace && trace.lgHeBug && !helike.lgHugeCaseB) )
01623 {
01624
01625
01626
01627
01628 fprintf(ioQQQ,"These indices may not represent the same levels for each element!\n");
01629 fprintf(ioQQQ,"index\tdesig\tE(wn,He)\tE(wn,O)\tE(wn,Fe)\n");
01630 for( ipHi=1; ipHi<( MAX3(iso.numLevels_max[ipHE_LIKE][ipHELIUM],
01631 iso.numLevels_max[ipHE_LIKE][ipIRON],iso.numLevels_max[ipHE_LIKE][ipOXYGEN]) - iso.nCollapsed_max[ipHE_LIKE][ipOXYGEN] ); ++ipHi )
01632 {
01633 float enHe=0. , enO=0., enFe=0.;
01634 nelem = -1;
01635 if( dense.lgElmtOn[ipHELIUM] &&(ipHi < iso.numLevels_max[ipHE_LIKE][ipHELIUM]) )
01636 {
01637 enHe = EmisLines[ipHE_LIKE][ipHELIUM][ipHi][0].EnergyWN;
01638 nelem = ipHELIUM;
01639 }
01640 else
01641 {
01642 enHe = 0.;
01643 }
01644 if( dense.lgElmtOn[ipOXYGEN] &&(ipHi < iso.numLevels_max[ipHE_LIKE][ipOXYGEN]) )
01645 {
01646 enO = EmisLines[ipHE_LIKE][ipOXYGEN][ipHi][0].EnergyWN;
01647 nelem = ipOXYGEN;
01648 }
01649 else
01650 {
01651 enO = 0.;
01652 }
01653 if( dense.lgElmtOn[ipIRON] &&(ipHi < iso.numLevels_max[ipHE_LIKE][ipIRON]) )
01654 {
01655 enFe = EmisLines[ipHE_LIKE][ipIRON][ipHi][0].EnergyWN;
01656 nelem = ipIRON;
01657 }
01658 else
01659 {
01660 enFe = 0.;
01661 }
01662 ASSERT( nelem >0 && nelem < LIMELM );
01663 if( ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
01664 {
01665 fprintf(ioQQQ,"%li\t%li %c%c\t%.5e\t%.5e\t%.5e\n",
01666 ipHi, N_(ipHi), chSpin[S_(ipHi)], chL[MIN2(5,L_(ipHi))],
01667 enHe, enO, enFe );
01668 }
01669
01670 }
01671
01672 fprintf(ioQQQ,"lo\thi\t He states \n");
01673 nelem = ipHELIUM;
01674 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem]; ++ipHi )
01675 {
01676 for( ipLo=ipHe1s1S; ipLo<ipHi; ++ipLo )
01677 {
01678 fprintf(ioQQQ," %li\t%li\t%li %c%c\t%li %c%c\n",
01679 ipLo, ipHi,
01680 N_(ipLo), chSpin[S_(ipLo)], chL[MIN2(5,L_(ipLo))],
01681 N_(ipHi), chSpin[S_(ipHi)], chL[MIN2(5,L_(ipHi))] );
01682 }
01683 }
01684 }
01685
01686 if( !helike.lgHugeCaseB )
01687 {
01688 {
01689
01690
01691 enum {DEBUG_LOC=false};
01692
01693 if( DEBUG_LOC )
01694 {
01695 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
01696 {
01697 if( dense.lgElmtOn[nelem] )
01698 {
01699 fprintf(ioQQQ,"%li\t%s\t",nelem+1,elementnames.chElementSym[nelem] );
01700
01701 ipLo = ipHe1s1S;
01702 for( ipHi=ipLo+1; ipHi<=ipHe2p1P; ++ipHi )
01703 {
01704 if( ipHi==ipHe2s1S || ipHi==ipHe2p3P0 ) continue;
01705 fprintf(ioQQQ,"%.4f\t", EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng);
01706 }
01707 fprintf(ioQQQ,"\n");
01708 }
01709 }
01710 cdEXIT(EXIT_SUCCESS);
01711 }
01712 }
01713
01714
01715
01716
01717
01718
01719 if( lgDataPathSet == true )
01720 {
01721
01722 strcpy( chFilename , chDataPath );
01723 strcat( chFilename , "he1_cs.dat" );
01724 }
01725 else
01726 {
01727
01728 strcpy( chFilename , "he1_cs.dat" );
01729 }
01730
01731 if( trace.lgTrace )
01732 fprintf( ioQQQ," HeCreate opening he1_cs.dat:");
01733
01734 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
01735 {
01736 fprintf( ioQQQ, " HeCreate could not open he1_cs.dat\n" );
01737 if( lgDataPathSet == true )
01738 fprintf( ioQQQ, " even tried path\n" );
01739
01740 if( lgDataPathSet == true )
01741 {
01742 fprintf( ioQQQ, " HeCreate could not open he1_cs.dat\n");
01743 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
01744 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
01745 }
01746
01747
01748 path_not_set();
01749
01750 puts( "[Stop in HeCreate]" );
01751 cdEXIT(EXIT_FAILURE);
01752 }
01753
01754
01755 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01756 {
01757 fprintf( ioQQQ, " HeCreate could not read first line of he1_cs.dat.\n");
01758 puts( "[Stop in HeCreate]" );
01759 cdEXIT(EXIT_FAILURE);
01760 }
01761 i = 1;
01762 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01763
01764
01765
01766
01767 if( i1 !=COLLISMAGIC )
01768 {
01769 fprintf( ioQQQ,
01770 " HeCreate: the version of he1_cs.dat is not the current version.\n" );
01771 fprintf( ioQQQ,
01772 " HeCreate: I expected to find the number %i and got %li instead.\n" ,
01773 COLLISMAGIC, i1 );
01774 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01775 puts( "[Stop in HeCreate]" );
01776 cdEXIT(EXIT_FAILURE);
01777 }
01778
01779
01780 lgHIT = false;
01781 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01782 {
01783
01784 if( chLine[0] != '#')
01785 {
01786 lgHIT = true;
01787 helike.nCS = 0;
01788 lgEOL = false;
01789 i = 1;
01790 while( !lgEOL && helike.nCS < HE1CSARRAY)
01791 {
01792 helike.CSTemp[helike.nCS] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01793 ++helike.nCS;
01794 }
01795 break;
01796 }
01797 }
01798 --helike.nCS;
01799 if( !lgHIT )
01800 {
01801 fprintf( ioQQQ, " HeCreate could not find line in CS temperatures in he1_cs.dat.\n");
01802 puts( "[Stop in HeCreate]" );
01803 cdEXIT(EXIT_FAILURE);
01804 }
01805
01806
01807 helike.HeCS = (float ****)MALLOC(sizeof(float ***)*(unsigned)LIMELM );
01808
01809 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
01810 {
01814 if( nelem != ipHELIUM )
01815 continue;
01816
01817
01818 if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
01819 {
01820 helike.HeCS[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]- iso.nCollapsed_max[ipHE_LIKE][nelem]) );
01821
01822
01823 helike.HeCS[nelem][0] = NULL;
01824 for( ipHi=ipHe2s3S; ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem];++ipHi )
01825 {
01826 helike.HeCS[nelem][ipHi] = (float**)MALLOC(sizeof(float*)*(unsigned)ipHi );
01827
01828 for( ipLo=ipHe1s1S; ipLo<ipHi; ++ipLo )
01829 {
01830 helike.HeCS[nelem][ipHi][ipLo] = (float*)MALLOC(sizeof(float)*(unsigned)helike.nCS );
01831
01832 for( i=0; i<helike.nCS; ++i )
01833 {
01834 helike.HeCS[nelem][ipHi][ipLo][i] = 0.f;
01835 }
01836 }
01837 }
01838 }
01839 }
01840
01841
01842 lgHIT = false;
01843 nelem = ipHELIUM;
01844 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01845 {
01846 char *chTemp;
01847
01848 if( (chLine[0] == ' ') || (chLine[0]=='\n') )
01849 break;
01850 if( chLine[0] != '#')
01851 {
01852 lgHIT = true;
01853
01854
01855 j = 1;
01856 ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
01857 ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
01858 ASSERT( ipLo < ipHi );
01859 if( ipHi >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
01860 continue;
01861 else
01862 {
01863 chTemp = chLine;
01864
01865 for( i=0; i<3; ++i )
01866 {
01867 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
01868 {
01869 fprintf( ioQQQ, " HeCreate no init cs\n" );
01870 puts( "[Stop in HeCreate]" );
01871 cdEXIT(EXIT_FAILURE);
01872 }
01873 ++chTemp;
01874 }
01875
01876 for( i=0; i<helike.nCS; ++i )
01877 {
01878 float a;
01879 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
01880 {
01881 fprintf( ioQQQ, " HeCreate not scan cs, current indices: %li %li\n", ipHi, ipLo );
01882 puts( "[Stop in HeCreate]" );
01883 cdEXIT(EXIT_FAILURE);
01884 }
01885 ++chTemp;
01886 sscanf( chTemp , "%e" , &a );
01887 helike.HeCS[nelem][ipHi][ipLo][i] = a;
01888 }
01889 }
01890 }
01891 }
01892
01893
01894 fclose( ioDATA );
01895
01896
01897
01898 {
01899
01900
01901
01902 enum {AGN=false};
01903
01904 if( AGN )
01905 {
01906 # define NTEMP 6
01907 double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
01908 double telog[NTEMP] ,
01909 cs ,
01910 ratecoef;
01911 nelem = ipHELIUM;
01912 fprintf( ioQQQ,"trans");
01913 for(i=0; i<NTEMP; ++i)
01914 {
01915 telog[i] = log10( te[i] );
01916 fprintf( ioQQQ,"\t%.3e",te[i]);
01917 }
01918 for(i=0; i<NTEMP; ++i)
01919 {
01920 fprintf( ioQQQ,"\t%.3e",te[i]);
01921 }
01922 fprintf(ioQQQ,"\n");
01923
01924 for( ipHi=ipHe2s3S; ipHi< iso.numLevels_max[ipHE_LIKE][ipHELIUM]; ++ipHi )
01925 {
01926 for( ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
01927 {
01928
01929 fprintf( ioQQQ,"%s - %s",
01930 chLevel[ipLo] , chLevel[ipHi] );
01931
01932
01933 for(i=0; i<NTEMP; ++i)
01934 {
01935 phycon.alogte = telog[i];
01936
01937 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
01938 fprintf(ioQQQ,"\t%.2e", cs );
01939 }
01940
01941
01942 for(i=0; i<NTEMP; ++i)
01943 {
01944 phycon.alogte = telog[i];
01945 phycon.te = pow(10.,telog[i] );
01946 tfidle(false);
01947 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
01948
01949 ratecoef = cs/sqrt(phycon.te)*COLL_CONST/iso.stat[ipHE_LIKE][nelem][ipLo] *
01950 sexp( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK / phycon.te );
01951 fprintf(ioQQQ,"\t%.2e", ratecoef );
01952 }
01953 fprintf(ioQQQ,"\n");
01954 }
01955 }
01956 cdEXIT(EXIT_FAILURE);
01957 }
01958 }
01959 }
01960
01961
01962
01963
01964
01965
01966 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
01967 {
01968 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
01969 {
01970 for( ipHi=ipHe2s3S; ipHi < iso.numLevels_max[ipHE_LIKE][nelem];++ipHi )
01971 {
01972 free( HeliumAs[nelem][ipHi] );
01973 }
01974 free( HeliumAs[nelem] );
01975 }
01976 }
01977 free( HeliumAs );
01978
01979 for( ipLo=ipHe1s1S; ipLo < max_num_levels;++ipLo )
01980 {
01981 free( RobbinsC[ipLo] );
01982 }
01983 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
01984 {
01985 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
01986 {
01987 free( energies[nelem] );
01988 free( SumAPerN[nelem] );
01989 }
01990 }
01991 free( energies );
01992 free( SumAPerN );
01993
01994 free( RobbinsC );
01995
01996
01997
01998
01999
02000
02001 if( helike.lgHugeCaseB )
02002 {
02003
02004 if( helike.lgRandErrGen )
02005 {
02006 HeLikeError(ipHELIUM);
02007 }
02008 HeRecom(ipHELIUM);
02009 }
02010
02011 DEBUG_EXIT( "HeCreate()" );
02012 return;
02013 }
02014