1 /***********************************************************************
2 
3  HiSIM (Hiroshima University STARC IGFET Model)
4  Copyright (C) 2014 Hiroshima University & STARC
5 
6  MODEL NAME : HiSIM_HV
7  ( VERSION : 2  SUBVERSION : 8  REVISION : 0 )
8  Model Parameter 'VERSION' : 2.80
9  FILE : hsm2eval_dep.h
10 
11  DATE : 2014.5.12
12 
13  released by
14                 Hiroshima University &
15                 Semiconductor Technology Academic Research Center (STARC)
16 ***********************************************************************/
17 
18 /**********************************************************************
19 
20 The following source code, and all copyrights, trade secrets or other
21 intellectual property rights in and to the source code in its entirety,
22 is owned by the Hiroshima University and the STARC organization.
23 
24 All users need to follow the "HISIM_HV Distribution Statement and
25 Copyright Notice" attached to HiSIM model.
26 
27 -----HISIM Distribution Statement and Copyright Notice--------------
28 
29 Software is distributed as is, completely without warranty or service
30 support. Hiroshima University or STARC and its employees are not liable
31 for the condition or performance of the software.
32 
33 Hiroshima University and STARC own the copyright and grant users a perpetual,
34 irrevocable, worldwide, non-exclusive, royalty-free license with respect
35 to the software as set forth below.
36 
37 Hiroshima University and STARC hereby disclaims all implied warranties.
38 
39 Hiroshima University and STARC grant the users the right to modify, copy,
40 and redistribute the software and documentation, both within the user's
41 organization and externally, subject to the following restrictions
42 
43 1. The users agree not to charge for Hiroshima University and STARC code
44 itself but may charge for additions, extensions, or support.
45 
46 2. In any product based on the software, the users agree to acknowledge
47 Hiroshima University and STARC that developed the software. This
48 acknowledgment shall appear in the product documentation.
49 
50 3. The users agree to reproduce any copyright notice which appears on
51 the software on any copy or modification of such made available
52 to others."
53 
54 Toshimasa Asahara, President, Hiroshima University
55 Mitiko Miura-Mattausch, Professor, Hiroshima University
56 Katsuhiro Shimohigashi, President&CEO, STARC
57 June 2008 (revised October 2011)
58 *************************************************************************/
59 
60 { // Begin : hsm2eval_dep
61 
62 /* define local variavles */
63   int    depmode ;
64   double afact, afact2, afact3, bfact, cfact ;
65   double W_bsub0, W_bsubL, W_s0, W_sL, W_sub0, W_subL, W_b0, W_bL, vthn ;
66   double phi_s0_DEP = 0.0, phi_sL_DEP = 0.0 , Vbi_DEP ;
67   double phi_s0_DEP_dVgs, phi_s0_DEP_dVbs, phi_s0_DEP_dVds ;
68   double phi_sL_DEP_dVgs, phi_sL_DEP_dVbs, phi_sL_DEP_dVds ;
69   double phi_j0_DEP, phi_jL_DEP, Psbmax, phi_b0_DEP_lim, phi_bL_DEP_lim ;
70 
71 
72   double phi_jL_DEP_dVgs, phi_jL_DEP_dVds, phi_jL_DEP_dVbs ;
73 
74   double Vgp0, Vgp1, Vgp0old, phi_j0_DEP_old, phi_jL_DEP_old, phi_b0_DEP_old, phi_bL_DEP_old, phi_s0_DEP_old, phi_sL_DEP_old ;
75   double phi_j0_DEP_acc, phi_jL_DEP_acc ;
76 
77 
78   double Q_s0, Q_sL = 0.0 ;
79   double Q_s0_dVgs, Q_sL_dVgs = 0.0, Q_s0_dVds, Q_sL_dVds = 0.0, Q_s0_dVbs, Q_sL_dVbs = 0.0 ;
80   double Q_sub0, Q_subL, Q_sub0_dVgs, Q_subL_dVgs, Q_sub0_dVds, Q_subL_dVds, Q_sub0_dVbs, Q_subL_dVbs ;
81   double Qn_res0, Qn_res0_dVgs, Qn_res0_dVds, Qn_res0_dVbs ;
82 
83 
84   double y1, y2, dety ;
85   double y11, y12 ;
86   double y21, y22 ;
87 
88   double y1_dVgs, y1_dVds, y1_dVbs ;
89   double y2_dVgs, y2_dVds, y2_dVbs ;
90 
91   double rev11 = 0.0, rev12 = 0.0 ;
92   double rev21 = 0.0, rev22 = 0.0 ;
93 
94   double phi_b0_DEP_ini ;
95   double y0, dydPsm ;
96 
97   double W_b0_dVgs, W_b0_dVds, W_b0_dVbs ;
98 
99   double W_res0 ;
100   double W_s0_dVgs, W_s0_dVds, W_s0_dVbs ;
101 
102   double phi_b0_DEP,  Q_b0_dep, Q_sub0_dep ;
103   double phi_b0_DEP_dVgs, phi_b0_DEP_dVds, phi_b0_DEP_dVbs ;
104   double phi_j0_DEP_dVgs, phi_j0_DEP_dVds, phi_j0_DEP_dVbs ;
105   double Q_b0_dep_dVgs, Q_b0_dep_dVds, Q_b0_dep_dVbs ;
106   double Q_sub0_dep_dVgs, Q_sub0_dep_dVds, Q_sub0_dep_dVbs ;
107 
108   double phi_bL_DEP, Q_bL_dep, Q_subL_dep ;
109   double phi_bL_DEP_dVgs, phi_bL_DEP_dVds, phi_bL_DEP_dVbs ;
110   double Q_bL_dep_dVgs, Q_bL_dep_dVds, Q_bL_dep_dVbs ;
111   double Q_subL_dep_dVgs, Q_subL_dep_dVds, Q_subL_dep_dVbs ;
112 
113   double q_Ndepm_esi, Idd_drift,Idd_diffu ;
114   double Qn_bac0 ;
115   double Qn_bac0_dVgs, Qn_bac0_dVds, Qn_bac0_dVbs ;
116 
117   double Mu_res, Mu_bac ;
118   double Mu_res_dVgs, Mu_res_dVds, Mu_res_dVbs ;
119   double Mu_bac_dVgs, Mu_bac_dVds, Mu_bac_dVbs ;
120 
121   double Q_n0_cur, Q_nL_cur ;
122   double Q_n0_cur_dVgs, Q_n0_cur_dVds, Q_n0_cur_dVbs ;
123   double Q_nL_cur_dVgs, Q_nL_cur_dVds, Q_nL_cur_dVbs ;
124 
125   double Q_s0_dep, Q_sL_dep ;
126   double Q_s0_dep_dVgs, Q_s0_dep_dVds, Q_s0_dep_dVbs ;
127   double Q_sL_dep_dVgs, Q_sL_dep_dVds, Q_sL_dep_dVbs ;
128 
129   double sm_delta ;
130   double phib_ref, phib_ref_dPs, phib_ref_dPd ;
131   double Q_s0_dPs, Q_sL_dPs, Q_s0_dPb, Q_sL_dPb ;
132   double Q_b0_dep_dPb, Q_bL_dep_dPb, Q_b0_dep_dPd, Q_bL_dep_dPd, Q_sub0_dep_dPd, Q_subL_dep_dPd ;
133   double phi_j0_DEP_dPb, phi_jL_DEP_dPb ;
134   double NdepmpNsub_inv1, NdepmpNsub ;
135 
136 
137 
138   double Q_n0, Q_n0_dVgs, Q_n0_dVds, Q_n0_dVbs ;
139   double Q_nL, Q_nL_dVgs, Q_nL_dVds, Q_nL_dVbs ;
140 
141   double phi_s0_DEP_ini, phi_sL_DEP_ini ;
142 
143 
144   double C_QE2, C_ESI2, Tn2 ;
145   double Ndepm2, q_Ndepm ;
146   double C_2ESIpq_Ndepm, C_2ESIpq_Ndepm_inv , C_2ESI_q_Ndepm ;
147   double C_2ESIpq_Nsub , C_2ESIpq_Nsub_inv ;
148   double ps_conv3 , ps_conv23 ;
149   double Ids_res, Ids_bac, Edri ;
150   double Ids_res_dVgs, Ids_res_dVds, Ids_res_dVbs ;
151   double Ids_bac_dVgs, Ids_bac_dVds, Ids_bac_dVbs ;
152   double Edri_dVgs, Edri_dVds, Edri_dVbs ;
153 
154   double T0_dVg, T0_dVb,T0_dVd ;
155   double T1_dVgs, T1_dVds, T1_dVbs ;
156   double T2_dVgs, T2_dVds, T2_dVbs ;
157   double T3_dVgs, T3_dVds, T3_dVbs ;
158   double T4_dVgs, T4_dVds, T4_dVbs ;
159   double T5_dVgs, T5_dVds, T5_dVbs ;
160 
161 
162   double Vgpp ;
163   double Vgpp_dVgs, Vgpp_dVds, Vgpp_dVbs ;
164   double Vdseff0, Vdseff0_dVgs, Vdseff0_dVds, Vdseff0_dVbs ;
165   double phib_ref_dVgs, phib_ref_dVds, phib_ref_dVbs ;
166 
167   double Qn_delta ;
168   double Qn_drift, Qn_drift_dVgs, Qn_drift_dVds, Qn_drift_dVbs ;
169 
170   double Ey_suf, Ey_suf_dVgs, Ey_suf_dVds, Ey_suf_dVbs ;
171 
172   double DEPQFN1 = 2.0 ;
173   double DEPQFN3 = 0.3 ;
174   double DEPQFN_dlt = 2.0 ;
175   double Ps_delta = 0.06 ;
176   double Ps_delta0 = 0.08 ;
177 
178 /*--------------------------------------*
179  * CeilingPow with derivatives for delta
180  *-----------------*/
181 #define Fn_SL_CP3( y , x , xmin , delta , pw , dx , dxmin, ddelta) { \
182  if(x < xmin + delta && delta >= 0.0) { \
183    double TMF0, TMF1; \
184    TMF1 = xmin + delta - x ; \
185    Fn_CP2( TMF0 , TMF1 , delta , pw , dx , ddelta )  \
186    y = xmin + delta - TMF0 ; \
187    dx = dx ; \
188    dxmin = 1.0-dx ; \
189    ddelta = 1.0-dx-ddelta; \
190  } else { \
191    y = x ; \
192    dx = 1.0 ; \
193    dxmin = 0.0 ; \
194    ddelta = 0.0 ; \
195  } \
196 }
197 
198 
199     // Constants
200     Vbi_DEP = here->HSM2_Vbipn ;
201     q_Ndepm = C_QE * here->HSM2_ndepm ;
202     Ndepm2  = here->HSM2_ndepm * here->HSM2_ndepm ;
203     q_Ndepm_esi = C_QE * here->HSM2_ndepm * C_ESI ;
204     q_Nsub = C_QE * here->HSM2_nsub ;
205     C_QE2  = C_QE * C_QE ;
206     C_ESI2 = C_ESI * C_ESI ;
207     Tn2    = model->HSM2_tndep * model->HSM2_tndep ;
208     C_2ESIpq_Ndepm = 2.0 * C_ESI/q_Ndepm ;
209     C_2ESIpq_Ndepm_inv = q_Ndepm / (2.0 * C_ESI) ;
210     C_2ESI_q_Ndepm = 2.0 * C_ESI * q_Ndepm ;
211     C_2ESIpq_Nsub  = 2.0 * C_ESI / q_Nsub  ;
212     C_2ESIpq_Nsub_inv  = q_Nsub / (2.0 * C_ESI) ;
213     NdepmpNsub  = here->HSM2_ndepm / here->HSM2_nsub ;
214     NdepmpNsub_inv1  = 1.0 / (1.0 + NdepmpNsub ) ;
215     ps_conv3  = ps_conv * 1000.0 ;
216     ps_conv23 = ps_conv2 * 1000.0 ;
217     here->HSM2_qnsub_esi = q_Ndepm_esi ;
218 
219 
220      //---------------------------------------------------*
221      // depletion MOS mode
222      //------------------//
223 
224      Vbsc = Vbs ;
225      Vbsc_dVbse = 1.0 ;
226 
227      /*---------------------------------------------------*
228       * initial potential phi_s0_DEP,phi_b0_DEP,phi_j0_DEP calculated.
229       *------------------*/
230 
231        Vgp = Vgp + epsm10 * 1.0e7 ;
232 
233 
234       afact = Cox * Cox / here->HSM2_cnst0 / here->HSM2_cnst0 ;
235       afact2 = afact / here->HSM2_nin / here->HSM2_nin * Ndepm2 ;
236       W_bsub0 = sqrt(2.0e0 * C_ESI / C_QE * here->HSM2_nsub / (here->HSM2_nsub
237               + here->HSM2_ndepm) / here->HSM2_ndepm * ( - Vbsc + Vbi_DEP)) ;
238 
239       if( W_bsub0 > model->HSM2_tndep ) {
240 
241         Vgp0 = 0.0;
242 
243         W_b0 = model->HSM2_tndep ;
244         phi_b0_DEP = 0.0 ;
245         phi_j0_DEP = phi_b0_DEP - C_2ESIpq_Ndepm_inv * W_b0 * W_b0 ;
246         phi_b0_DEP_lim = 0.0 ;
247 
248         Vgp0old = Vgp0 ;
249         phi_j0_DEP_old = phi_j0_DEP ;
250 
251         for ( lp_s0 = 1 ; lp_s0 <= lp_se_max + 1 ; lp_s0 ++ ) {
252 
253           W_b0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_j0_DEP) ) ;
254           Fn_SU_CP( W_b0 , W_b0 , model->HSM2_tndep , 1e-8, 2 , T0 )
255           W_sub0 = sqrt(C_2ESIpq_Nsub * (phi_j0_DEP - Vbsc + Vbi_DEP) ) ;
256 
257           Q_b0_dep = W_b0 * q_Ndepm ;
258           Q_b0_dep_dPd = - C_ESI / W_b0 * T0 ;
259           Q_sub0_dep = - W_sub0 * q_Nsub ;
260           Q_sub0_dep_dPd = - C_ESI / W_sub0 ;
261 
262           y1 = Cox * (Vgp0 - phi_b0_DEP) + Q_b0_dep + Q_sub0_dep ;
263           y11 = Cox ;
264           y12 = Q_b0_dep_dPd + Q_sub0_dep_dPd ;
265 
266           y2 = phi_j0_DEP - NdepmpNsub_inv1 * (NdepmpNsub * phi_b0_DEP + Vbsc - Vbi_DEP) ;
267           y21 = 0.0 ;
268           y22 = 1.0 ;
269 
270           dety = y11 * y22 - y21 * y12;
271           rev11 = (y22) / dety ;
272           rev12 = ( - y12) / dety ;
273           rev21 = ( - y21) / dety ;
274           rev22 = (y11) / dety ;
275 
276           if( fabs( rev11 * y1 + rev12 * y2 ) > 0.5 ) {
277             Vgp0 = Vgp0 - 0.5 * Fn_Sgn( rev11 * y1 + rev12 * y2 ) ;
278             phi_j0_DEP = phi_j0_DEP - 0.5 * Fn_Sgn( rev21 * y1 + rev22 * y2 ) ;
279           } else {
280             Vgp0 = Vgp0 - ( rev11 * y1 + rev12 * y2 ) ;
281             phi_j0_DEP = phi_j0_DEP - ( rev21 * y1 + rev22 * y2 ) ;
282           }
283 
284           if( fabs(Vgp0 - Vgp0old) <= ps_conv &&
285                 fabs(phi_j0_DEP - phi_j0_DEP_old) <= ps_conv ) lp_s0=lp_se_max + 1 ;
286 
287           Vgp0old = Vgp0 ;
288           phi_j0_DEP_old = phi_j0_DEP ;
289         }
290         phi_j0_DEP_acc = phi_j0_DEP ;
291 
292         W_sub0 = model->HSM2_tndep * NdepmpNsub ;
293         phi_j0_DEP = C_2ESIpq_Nsub_inv * W_sub0 * W_sub0 + Vbsc - Vbi_DEP ;
294         phi_b0_DEP = phi_j0_DEP + C_2ESIpq_Ndepm_inv * Tn2 ;
295         phi_s0_DEP = phi_b0_DEP ;
296         Psbmax = phi_b0_DEP ;
297         Vgp1 = phi_b0_DEP ;
298         if( Vgp > Vgp0 ) {
299           depmode = 1 ;
300         } else if(Vgp > Vgp1 ) {
301           depmode = 3 ;
302         } else {
303           depmode = 2 ;
304         }
305 
306       } else {
307         Vgp0 = 0.0 ;
308         Vgp1 = Vgp0 ;
309         Psbmax = 0.0 ;
310         phi_b0_DEP_lim = Vgp0 ;
311         W_b0 = W_bsub0 ;
312         W_sub0 = W_b0 * NdepmpNsub ;
313         phi_j0_DEP = C_2ESIpq_Nsub_inv * W_sub0 * W_sub0 + Vbsc - Vbi_DEP ;
314         phi_b0_DEP = C_2ESIpq_Ndepm_inv * W_b0 * W_b0 + phi_j0_DEP ;
315         phi_j0_DEP_acc = phi_j0_DEP ;
316         if( Vgp > Vgp0 ) {
317           depmode = 1 ;
318         } else {
319           depmode = 2 ;
320         }
321 
322       }
323 
324 
325       T1 = C_2ESI_q_Ndepm * ( Psbmax - ( - here->HSM2_Pb2n + Vbsc)) ;
326       if ( T1 > 0.0 ) {
327         vthn = - here->HSM2_Pb2n + Vbsc - sqrt(T1) / Cox ;
328       } else {
329         vthn = - here->HSM2_Pb2n + Vbsc ;
330       }
331 
332       /* primary value */
333 
334        if( Vgp > Vgp0 ) {
335         /* accumulation region */
336          phi_j0_DEP = phi_j0_DEP_acc ;
337          phi_b0_DEP = 0.0 ;
338          phi_s0_DEP_ini = log(afact * Vgp * Vgp) / (beta + 2.0 / Vgp) + phi_b0_DEP ;
339 
340          if( phi_s0_DEP_ini < phi_b0_DEP_lim + ps_conv23 ) phi_s0_DEP_ini = phi_b0_DEP_lim + ps_conv23 ;
341 
342        } else if( Vgp > Vgp1 ) {
343         /* depletion region */
344 
345          phi_s0_DEP_ini = phi_s0_DEP ;
346 
347        } else {
348         /* depletion and inversion region */
349 
350          if( Vgp > vthn ) {
351           /* depletion */
352            bfact = - 2.0 * afact * Vgp + beta ;
353            cfact = afact * Vgp * Vgp - beta * phi_b0_DEP ;
354            phi_b0_DEP_old = phi_b0_DEP ;
355 
356            phi_s0_DEP_ini = ( - bfact + sqrt(bfact * bfact - 4.0 * afact * cfact)) / 2.0 / afact ;
357            if( phi_s0_DEP_ini > Psbmax - ps_conv3 ) phi_s0_DEP_ini = Psbmax - ps_conv3 ;
358 
359            W_s0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_s0_DEP_ini) ) ;
360            W_b0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_j0_DEP) ) ;
361 
362            if( W_s0 + W_b0 > model->HSM2_tndep ) {
363              for ( lp_s0 = 1 ; lp_s0 <= lp_se_max + 1 ; lp_s0 ++ ) {
364 
365                y0 = W_s0 + W_b0 - model->HSM2_tndep ;
366 
367                dydPsm = C_ESI / q_Ndepm / W_s0
368                   + C_ESI / q_Ndepm * ( 1.0 - (here->HSM2_ndepm
369                   / here->HSM2_nsub) / ( 1.0 + (NdepmpNsub))) / W_b0 ;
370 
371                if( fabs(y0 / dydPsm) > 0.5 ) {
372                  phi_b0_DEP = phi_b0_DEP - 0.5 * Fn_Sgn(y0 / dydPsm) ;
373                } else {
374                  phi_b0_DEP = phi_b0_DEP - y0 / dydPsm ;
375                }
376 
377                if( (phi_b0_DEP - Vbsc + Vbi_DEP) < epsm10 )
378                   phi_b0_DEP=Vbsc - Vbi_DEP + epsm10 ;
379 
380                cfact = afact * Vgp * Vgp - beta * phi_b0_DEP ;
381                T1 = bfact * bfact - 4.0 * afact * cfact ;
382                if( T1 > 0.0 ) {
383                  phi_s0_DEP_ini = ( - bfact + sqrt(T1)) / 2.0 / afact ;
384                } else {
385                  phi_s0_DEP_ini = ( - bfact) / 2.0 / afact ;
386                }
387 
388                if( phi_s0_DEP_ini > Psbmax ) phi_s0_DEP_ini = Psbmax ;
389                if( phi_s0_DEP_ini > phi_b0_DEP ) {
390                  phi_s0_DEP_ini = phi_b0_DEP - ps_conv23 ;
391                  lp_s0=lp_se_max + 1 ;
392                }
393 
394                W_s0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_s0_DEP_ini) ) ;
395                phi_j0_DEP = ( NdepmpNsub * phi_b0_DEP
396                  + Vbsc - Vbi_DEP) / (1.0 + NdepmpNsub) ;
397                W_b0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_j0_DEP) ) ;
398 
399               if( fabs(phi_b0_DEP - phi_b0_DEP_old) <= 1.0e-8 ) lp_s0=lp_se_max + 1 ;
400                phi_b0_DEP_old = phi_b0_DEP ;
401              }
402            }
403 
404          } else {
405            afact3 = afact2 / exp(beta * Vbsc) ;
406            phi_b0_DEP_old = phi_b0_DEP ;
407            phi_s0_DEP_ini = log(afact3 * Vgp * Vgp) / ( - beta + 2.0 / Vgp) ;
408            W_s0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_s0_DEP_ini) ) ;
409            W_b0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_j0_DEP) ) ;
410            if( W_s0 + W_b0 >  model->HSM2_tndep ) {
411              for ( lp_s0 = 1 ; lp_s0 <= lp_s0_max + 1 ; lp_s0 ++ ) {
412 
413                y0 = W_s0 + W_b0 - model->HSM2_tndep ;
414                dydPsm = C_ESI / q_Ndepm / W_s0
415                  + C_ESI / q_Ndepm * ( 1.0 - (here->HSM2_ndepm /
416                  here->HSM2_nsub) / ( 1.0 + (NdepmpNsub))) / W_b0 ;
417 
418                if( fabs(y0 / dydPsm) > 0.5 ) {
419                  phi_b0_DEP = phi_b0_DEP - 0.5 * Fn_Sgn(y0 / dydPsm) ;
420                } else {
421                  phi_b0_DEP = phi_b0_DEP - y0 / dydPsm ;
422                }
423                if( (phi_b0_DEP - Vbsc + Vbi_DEP) < epsm10 )
424                     phi_b0_DEP=Vbsc - Vbi_DEP + epsm10 ;
425 
426                W_s0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_s0_DEP_ini) ) ;
427                phi_j0_DEP = ( NdepmpNsub * phi_b0_DEP
428                  + Vbsc - Vbi_DEP) / (1.0 + NdepmpNsub) ;
429                W_b0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_j0_DEP) ) ;
430 
431                if( fabs(phi_b0_DEP - phi_b0_DEP_old) <= 1.0e-5 ) lp_s0=lp_s0_max + 1 ;
432                phi_b0_DEP_old = phi_b0_DEP ;
433              }
434 
435            }
436          } // end of phi_b0_DEP loop //
437 
438        }
439        phi_b0_DEP_ini = phi_b0_DEP ;
440 
441        /*                               */
442        /* solve poisson at source side  */
443        /*                               */
444 
445        sm_delta = 0.12 ;
446 
447        flg_conv = 0 ;
448 
449        phi_s0_DEP = phi_s0_DEP_ini ;
450        phi_b0_DEP = phi_b0_DEP_ini ;
451 
452        phi_s0_DEP_old = phi_s0_DEP ;
453        phi_b0_DEP_old = phi_b0_DEP ;
454 
455        for ( lp_s0 = 1 ; lp_s0 <= lp_se_max + 1 ; lp_s0 ++ ) {
456 
457          phi_j0_DEP  = NdepmpNsub_inv1 * (NdepmpNsub * phi_b0_DEP + Vbsc - Vbi_DEP) ;
458          phi_j0_DEP_dPb = NdepmpNsub_inv1 * NdepmpNsub ;
459 
460          T1 = phi_b0_DEP - phi_j0_DEP ;
461          Fn_SL_CP( T2 , T1 , 0.0 , 0.1, 2 , T7 )
462          W_b0 = sqrt(C_2ESIpq_Ndepm * (T2) ) ;
463          Fn_SU_CP( W_b0 , W_b0 , model->HSM2_tndep, 1e-8, 2 , T8 )
464          W_sub0 = sqrt(C_2ESIpq_Nsub * (phi_j0_DEP - Vbsc + Vbi_DEP) ) ;
465          Q_b0_dep = W_b0 * q_Ndepm ;
466          Q_b0_dep_dPb = C_ESI / W_b0 * T7 * T8 ;
467          Q_b0_dep_dPd = - C_ESI / W_b0 * T7 * T8 ;
468          Q_sub0_dep = - W_sub0 * q_Nsub ;
469          Q_sub0_dep_dPd = - C_ESI / W_sub0 ;
470 
471          T10 = 8.0 * q_Ndepm_esi * Tn2 ;
472          phib_ref = (4.0 * phi_j0_DEP * phi_j0_DEP * C_ESI2 - 8.0 * phi_j0_DEP * C_ESI2 * phi_s0_DEP
473                   + 4.0 * C_ESI2 * phi_s0_DEP * phi_s0_DEP
474                   + 4.0 * phi_j0_DEP * q_Ndepm_esi * Tn2
475                   + 4.0 * phi_s0_DEP * q_Ndepm_esi * Tn2
476                   + Ndepm2 * C_QE2 * Tn2 * Tn2 ) / T10 ;
477          phib_ref_dPs = ( - 8.0 * phi_j0_DEP * C_ESI2 + 8.0 * C_ESI2 * phi_s0_DEP
478                       + 4.0 * q_Ndepm_esi * Tn2 ) / T10 ;
479          phib_ref_dPd = (   8.0 * phi_j0_DEP * C_ESI2 - 8.0 * C_ESI2 * phi_s0_DEP
480                       + 4.0 * q_Ndepm_esi * Tn2 ) / T10 ;
481 
482          T1 = beta * (phi_s0_DEP - phi_b0_DEP) ;
483          T2 = exp(T1) ;
484          if( phi_s0_DEP >= phi_b0_DEP ) {
485            Q_s0 = - here->HSM2_cnst0 * sqrt(T2 - 1.0 - T1 + 1e-15) ;
486            Q_s0_dPs = 0.5 * here->HSM2_cnst0 * here->HSM2_cnst0 / Q_s0 * (beta * T2 - beta ) ;
487            Q_s0_dPb = - Q_s0_dPs ;
488          } else {
489            T3 = exp( - beta * (phi_s0_DEP - Vbsc)) ;
490            T4 = exp( - beta * (phi_b0_DEP - Vbsc)) ;
491            Q_s0 = here->HSM2_cnst0 * sqrt(T2 - 1.0 - T1 + 1e-15 + here->HSM2_cnst1 * (T3 - T4) ) ;
492            T5 = 0.5 * here->HSM2_cnst0 * here->HSM2_cnst0 / Q_s0 ;
493            Q_s0_dPs = T5 * (beta * T2 - beta + here->HSM2_cnst1 * ( - beta * T3) ) ;
494            Q_s0_dPb = T5 * ( - beta * T2 + beta + here->HSM2_cnst1 * beta * T4 ) ;
495          }
496 
497          Fn_SU_CP( T1 , phib_ref , phi_b0_DEP_lim , sm_delta, 4 , T9 )
498 
499          y1 = phi_b0_DEP - T1 ;
500          y11 = - phib_ref_dPs * T9 ;
501          y12 = 1.0 - phib_ref_dPd * phi_j0_DEP_dPb * T9 ;
502 
503          y2 = Cox * (Vgp - phi_s0_DEP) + Q_s0 + Q_b0_dep + Q_sub0_dep ;
504          y21 = - Cox + Q_s0_dPs ;
505          y22 = Q_s0_dPb + Q_b0_dep_dPb + Q_b0_dep_dPd * phi_j0_DEP_dPb + Q_sub0_dep_dPd * phi_j0_DEP_dPb ;
506 
507          dety = y11 * y22 - y21 * y12;
508          rev11 = (y22) / dety ;
509          rev12 = ( - y12) / dety ;
510          rev21 = ( - y21) / dety ;
511          rev22 = (y11) / dety ;
512          if( fabs( rev21 * y1 + rev22 * y2 ) > 0.5 ) {
513            phi_s0_DEP = phi_s0_DEP - 0.5 * Fn_Sgn( rev11 * y1 + rev12 * y2 ) ;
514            phi_b0_DEP = phi_b0_DEP - 0.5 * Fn_Sgn( rev21 * y1 + rev22 * y2 ) ;
515          } else {
516            phi_s0_DEP = phi_s0_DEP - ( rev11 * y1 + rev12 * y2 ) ;
517            phi_b0_DEP = phi_b0_DEP - ( rev21 * y1 + rev22 * y2 ) ;
518          }
519 
520          if( fabs(phi_s0_DEP - phi_s0_DEP_old) <= ps_conv &&  fabs(phi_b0_DEP - phi_b0_DEP_old) <= ps_conv ) {
521            lp_s0=lp_se_max + 1 ;
522            flg_conv = 1 ;
523          }
524 
525          phi_s0_DEP_old = phi_s0_DEP ;
526          phi_b0_DEP_old = phi_b0_DEP ;
527 
528        }
529 
530        if( flg_conv == 0 ) {
531          printf( "*** warning(HiSIM(%s)): Went Over Iteration Maximum(Ps0)\n",model->HSM2modName ) ;
532          printf( " Vbse   = %7.3f Vdse = %7.3f Vgse = %7.3f\n" ,Vbse , Vdse , Vgse ) ;
533        }
534 
535       /* caluculate derivative */
536 
537        y1_dVgs = 0.0 ;
538        y1_dVds = 0.0 ;
539        y1_dVbs = - (8.0 * phi_j0_DEP * C_ESI2 - 8.0 * C_ESI2 * phi_s0_DEP
540                   + 4.0 * q_Ndepm_esi * Tn2) / T10
541                   * T9 * NdepmpNsub_inv1 * Vbsc_dVbse ;
542 
543        Q_b0_dep_dVbs = - C_ESI / W_b0 * T7 * T8 * NdepmpNsub_inv1 * Vbsc_dVbse ;
544 
545        Q_sub0_dep_dVbs = - C_ESI / W_sub0 * (NdepmpNsub_inv1 * Vbsc_dVbse - Vbsc_dVbse) ;
546 
547        T1 = beta * (phi_s0_DEP - phi_b0_DEP) ;
548        T2 = exp(T1) ;
549        if( phi_s0_DEP >= phi_b0_DEP ) {
550          Q_s0_dVbs = 0.0 ;
551        } else {
552          T3 = exp( - beta * (phi_s0_DEP - Vbsc)) ;
553          T4 = exp( - beta * (phi_b0_DEP - Vbsc)) ;
554          T5 = sqrt(T2 - 1.0 - T1 + 1e-15 + here->HSM2_cnst1 * (T3 - T4)) ;
555          Q_s0_dVbs = here->HSM2_cnst0 / 2.0 / T5 *
556                   (here->HSM2_cnst1 * (beta * T3 * Vbsc_dVbse - beta * T4 * Vbsc_dVbse) ) ;
557        }
558 
559        y2_dVgs = Cox_dVg * (Vgp - phi_s0_DEP) + Cox * Vgp_dVgs ;
560        y2_dVds = Cox_dVd * (Vgp - phi_s0_DEP) + Cox * Vgp_dVds ;
561        y2_dVbs = Cox_dVb * (Vgp - phi_s0_DEP) + Cox * Vgp_dVbs + Q_s0_dVbs + Q_b0_dep_dVbs + Q_sub0_dep_dVbs ;
562 
563        phi_s0_DEP_dVgs = - ( rev11 * y1_dVgs + rev12 * y2_dVgs ) ;
564        phi_s0_DEP_dVds = - ( rev11 * y1_dVds + rev12 * y2_dVds ) ;
565        phi_s0_DEP_dVbs = - ( rev11 * y1_dVbs + rev12 * y2_dVbs ) ;
566 
567        phi_b0_DEP_dVgs = - ( rev21 * y1_dVgs + rev22 * y2_dVgs ) ;
568        phi_b0_DEP_dVds = - ( rev21 * y1_dVds + rev22 * y2_dVds ) ;
569        phi_b0_DEP_dVbs = - ( rev21 * y1_dVbs + rev22 * y2_dVbs ) ;
570 
571        if( W_bsub0 > model->HSM2_tndep && depmode !=2 ) {
572          Fn_SU_CP2(phi_b0_DEP , phi_b0_DEP , phi_s0_DEP , 0.04, 2 , T1, T2 )  // HV_dlt=0.02
573          phi_b0_DEP_dVgs = phi_b0_DEP_dVgs * T1 + phi_s0_DEP_dVgs * T2 ;
574          phi_b0_DEP_dVds = phi_b0_DEP_dVds * T1 + phi_s0_DEP_dVds * T2 ;
575          phi_b0_DEP_dVbs = phi_b0_DEP_dVbs * T1 + phi_s0_DEP_dVbs * T2 ;
576        }
577 
578        phi_j0_DEP  = NdepmpNsub_inv1 * (NdepmpNsub * phi_b0_DEP + Vbsc - Vbi_DEP) ;
579        phi_j0_DEP_dVgs = NdepmpNsub_inv1 * NdepmpNsub * phi_b0_DEP_dVgs ;
580        phi_j0_DEP_dVds = NdepmpNsub_inv1 * NdepmpNsub * phi_b0_DEP_dVds ;
581        phi_j0_DEP_dVbs = NdepmpNsub_inv1 * NdepmpNsub * phi_b0_DEP_dVbs + NdepmpNsub_inv1 * Vbsc_dVbse  ;
582 
583        phib_ref = (4.0 * phi_j0_DEP * phi_j0_DEP * C_ESI2 - 8.0 * phi_j0_DEP * C_ESI2 * phi_s0_DEP
584                 + 4.0 * C_ESI2 * phi_s0_DEP * phi_s0_DEP
585                 + 4.0 * phi_j0_DEP * q_Ndepm_esi * Tn2
586                 + 4.0 * phi_s0_DEP * q_Ndepm_esi * Tn2
587                 + Ndepm2 * C_QE2 * Tn2 * Tn2 ) / T10 ;
588 
589        phib_ref_dVgs = ( 8.0 * phi_j0_DEP * phi_j0_DEP_dVgs * C_ESI2 - 8.0 * phi_j0_DEP_dVgs * C_ESI2 * phi_s0_DEP
590                        - 8.0 * phi_j0_DEP * C_ESI2 * phi_s0_DEP_dVgs + 8.0 * C_ESI2 * phi_s0_DEP * phi_s0_DEP_dVgs
591                      + 4.0 * phi_j0_DEP_dVgs * q_Ndepm_esi * Tn2
592                      + 4.0 * phi_s0_DEP_dVgs * q_Ndepm_esi * Tn2 ) / T10 ;
593        phib_ref_dVds = ( 8.0 * phi_j0_DEP * phi_j0_DEP_dVds * C_ESI2 - 8.0 * phi_j0_DEP_dVds * C_ESI2 * phi_s0_DEP
594                        - 8.0 * phi_j0_DEP * C_ESI2 * phi_s0_DEP_dVds + 8.0 * C_ESI2 * phi_s0_DEP * phi_s0_DEP_dVds
595                      + 4.0 * phi_j0_DEP_dVds * q_Ndepm_esi * Tn2
596                      + 4.0 * phi_s0_DEP_dVds * q_Ndepm_esi * Tn2 ) / T10 ;
597        phib_ref_dVbs = ( 8.0 * phi_j0_DEP * phi_j0_DEP_dVbs * C_ESI2 - 8.0 * phi_j0_DEP_dVbs * C_ESI2 * phi_s0_DEP
598                        - 8.0 * phi_j0_DEP * C_ESI2 * phi_s0_DEP_dVbs + 8.0 * C_ESI2 * phi_s0_DEP * phi_s0_DEP_dVbs
599                      + 4.0 * phi_j0_DEP_dVbs * q_Ndepm_esi * Tn2
600                      + 4.0 * phi_s0_DEP_dVbs * q_Ndepm_esi * Tn2 ) / T10 ;
601 
602        T1 = beta * (phi_s0_DEP - phi_b0_DEP) ;
603        T1_dVgs = beta * (phi_s0_DEP_dVgs - phi_b0_DEP_dVgs) ;
604        T1_dVds = beta * (phi_s0_DEP_dVds - phi_b0_DEP_dVds) ;
605        T1_dVbs = beta * (phi_s0_DEP_dVbs - phi_b0_DEP_dVbs) ;
606 
607        T2 = exp(T1) ;
608        T2_dVgs = T1_dVgs * T2 ;
609        T2_dVds = T1_dVds * T2 ;
610        T2_dVbs = T1_dVbs * T2 ;
611 
612        if( phi_s0_DEP >= phi_b0_DEP ) {
613 
614          T3 = sqrt(T2 - 1.0e0 - T1 + 1e-15 ) ;
615          T3_dVgs = (T2_dVgs - T1_dVgs) / 2.0 / T3 ;
616          T3_dVds = (T2_dVds - T1_dVds) / 2.0 / T3 ;
617          T3_dVbs = (T2_dVbs - T1_dVbs) / 2.0 / T3 ;
618 
619          Q_s0 = - here->HSM2_cnst0 * T3 ;
620 
621          Q_s0_dep = 0.0 ;
622          Q_sub0 = 0.0 ;
623 
624          W_b0 = sqrt(C_2ESIpq_Ndepm * (phi_b0_DEP - phi_j0_DEP) ) ;
625          Fn_SU_CP( T9 , W_b0 , model->HSM2_tndep, 5e-8, 2 , T4 )  // HV_dlt=1e-8
626 
627          W_sub0 = sqrt(C_2ESIpq_Nsub * (phi_j0_DEP - Vbsc + Vbi_DEP) ) ;
628          Q_b0_dep = T9 * q_Ndepm ;
629          Q_sub0_dep = - W_sub0 * q_Nsub ;
630 
631         /* derivative */
632          Q_s0_dVgs = - here->HSM2_cnst0 * T3_dVgs ;
633          Q_s0_dVds = - here->HSM2_cnst0 * T3_dVds ;
634          Q_s0_dVbs = - here->HSM2_cnst0 * T3_dVbs ;
635 
636          Q_n0 = Q_s0 ;
637          Q_n0_dVgs = Q_s0_dVgs ;
638          Q_n0_dVds = Q_s0_dVds ;
639          Q_n0_dVbs = Q_s0_dVbs ;
640 
641          Q_b0_dep_dVgs = C_ESI / W_b0 * (phi_b0_DEP_dVgs - phi_j0_DEP_dVgs) * T4 ;
642          Q_b0_dep_dVds = C_ESI / W_b0 * (phi_b0_DEP_dVds - phi_j0_DEP_dVds) * T4 ;
643          Q_b0_dep_dVbs = C_ESI / W_b0 * (phi_b0_DEP_dVbs - phi_j0_DEP_dVbs) * T4 ;
644 
645          Q_sub0_dep_dVgs = - C_ESI / W_sub0 * phi_j0_DEP_dVgs ;
646          Q_sub0_dep_dVds = - C_ESI / W_sub0 * phi_j0_DEP_dVds ;
647          Q_sub0_dep_dVbs = - C_ESI / W_sub0 * (phi_j0_DEP_dVbs - Vbsc_dVbse) ;
648 
649          Q_sub0_dVgs = 0.0 ;
650          Q_sub0_dVds = 0.0 ;
651          Q_sub0_dVbs = 0.0 ;
652 
653          Q_s0_dep_dVgs = 0.0 ;
654          Q_s0_dep_dVds = 0.0 ;
655          Q_s0_dep_dVbs = 0.0 ;
656 
657        } else {
658 
659          T3 = exp( - beta * (phi_s0_DEP - Vbsc)) ;
660          T4 = exp( - beta * (phi_b0_DEP - Vbsc)) ;
661          T5 = sqrt(T2 - 1.0 - T1 + here->HSM2_cnst1 * (T3 - T4) + 1e-15) ;
662          Q_s0 = here->HSM2_cnst0 * T5 ;
663 
664          T3_dVgs = - beta * T3 * phi_s0_DEP_dVgs ;
665          T3_dVds = - beta * T3 * phi_s0_DEP_dVds ;
666          T3_dVbs = - beta * T3 * (phi_s0_DEP_dVbs - Vbsc_dVbse) ;
667 
668          T4_dVgs = - beta * T4 * phi_b0_DEP_dVgs ;
669          T4_dVds = - beta * T4 * phi_b0_DEP_dVds ;
670          T4_dVbs = - beta * T4 * (phi_b0_DEP_dVbs - Vbsc_dVbse) ;
671 
672          T5_dVgs = (T2_dVgs - T1_dVgs + here->HSM2_cnst1 * (T3_dVgs - T4_dVgs)) / 2.0 / T5 ;
673          T5_dVds = (T2_dVds - T1_dVds + here->HSM2_cnst1 * (T3_dVds - T4_dVds)) / 2.0 / T5 ;
674          T5_dVbs = (T2_dVbs - T1_dVbs + here->HSM2_cnst1 * (T3_dVbs - T4_dVbs)) / 2.0 / T5 ;
675 
676          Q_s0_dVgs = here->HSM2_cnst0 * T5_dVgs ;
677          Q_s0_dVds = here->HSM2_cnst0 * T5_dVds ;
678          Q_s0_dVbs = here->HSM2_cnst0 * T5_dVbs ;
679 
680          if( W_bsub0 > model->HSM2_tndep && depmode !=2 ) {
681            Q_sub0 = 0.0 ;
682            Q_s0_dep = 0.0 ;
683 
684            Q_sub0_dVgs = 0.0 ;
685            Q_sub0_dVds = 0.0 ;
686            Q_sub0_dVbs = 0.0 ;
687 
688            Q_s0_dep_dVgs = 0.0 ;
689            Q_s0_dep_dVds = 0.0 ;
690            Q_s0_dep_dVbs = 0.0 ;
691          } else {
692            T3 = exp( - beta * (phi_s0_DEP - Vbsc)) ;
693            T4 = exp( - beta * (phi_b0_DEP - Vbsc)) ;
694            T5 = sqrt( - T1 + here->HSM2_cnst1 * (T3 - T4)) ;
695            Q_sub0 = here->HSM2_cnst0 * T5 - here->HSM2_cnst0 * sqrt( - T1)  ;
696            T6 = sqrt(T2 - 1.0e0 - T1 + 1e-15) ;
697            Q_s0_dep = here->HSM2_cnst0 * T6 ;
698 
699            Q_sub0_dVgs = here->HSM2_cnst0 / 2.0 / T5 * ( - beta * (phi_s0_DEP_dVgs - phi_b0_DEP_dVgs)
700               + here->HSM2_cnst1 * ( - beta * T3 * phi_s0_DEP_dVgs + beta * T4 * phi_b0_DEP_dVgs))
701               - here->HSM2_cnst0 / 2.0 / sqrt( - T1) * ( - beta * (phi_s0_DEP_dVgs - phi_b0_DEP_dVgs)) ;
702            Q_sub0_dVds = here->HSM2_cnst0 / 2.0 / T5 * ( - beta * (phi_s0_DEP_dVds - phi_b0_DEP_dVds)
703               + here->HSM2_cnst1 * ( - beta * T3 * phi_s0_DEP_dVds + beta * T4 * phi_b0_DEP_dVds))
704               - here->HSM2_cnst0 / 2.0 / sqrt( - T1) * ( - beta * (phi_s0_DEP_dVds - phi_b0_DEP_dVds)) ;
705            Q_sub0_dVbs = here->HSM2_cnst0 / 2.0 / T5 * ( - beta * (phi_s0_DEP_dVbs - phi_b0_DEP_dVbs)
706               + here->HSM2_cnst1 * ( - beta * T3 * (phi_s0_DEP_dVbs - Vbsc_dVbse) + beta * T4 * (phi_b0_DEP_dVbs - Vbsc_dVbse)))
707               - here->HSM2_cnst0 / 2.0 / sqrt( - T1) * ( - beta * (phi_s0_DEP_dVbs - phi_b0_DEP_dVbs)) ;
708 
709            Q_s0_dep_dVgs = here->HSM2_cnst0 / 2.0 / T6 * beta * (phi_s0_DEP_dVgs - phi_b0_DEP_dVgs) * (T2 - 1) ;
710            Q_s0_dep_dVds = here->HSM2_cnst0 / 2.0 / T6 * beta * (phi_s0_DEP_dVds - phi_b0_DEP_dVds) * (T2 - 1) ;
711            Q_s0_dep_dVbs = here->HSM2_cnst0 / 2.0 / T6 * beta * (phi_s0_DEP_dVbs - phi_b0_DEP_dVbs) * (T2 - 1) ;
712 
713          }
714 
715          Q_n0 = 0.0 ;
716          Q_n0_dVgs = 0.0 ;
717          Q_n0_dVds = 0.0 ;
718          Q_n0_dVbs = 0.0 ;
719 
720 
721          T1 = phi_b0_DEP - phi_j0_DEP ;
722          Fn_SL_CP( T2 , T1 , 0.0 , 0.1, 2 , T0 )
723          W_b0 = sqrt(C_2ESIpq_Ndepm * (T2) ) ;
724          Fn_SU_CP( T9 , W_b0 , model->HSM2_tndep, 5e-8, 2 , T3 )  // HV_dlt=1e-8
725          W_sub0 = sqrt(C_2ESIpq_Nsub * (phi_j0_DEP - Vbsc + Vbi_DEP) ) ;
726          Q_b0_dep = T9 * q_Ndepm ;
727          Q_sub0_dep = - W_sub0 * q_Nsub ;
728 
729          Q_b0_dep_dVgs = C_ESI / W_b0 * (phi_b0_DEP_dVgs - phi_j0_DEP_dVgs) * T0 * T3 ;
730          Q_b0_dep_dVds = C_ESI / W_b0 * (phi_b0_DEP_dVds - phi_j0_DEP_dVds) * T0 * T3 ;
731          Q_b0_dep_dVbs = C_ESI / W_b0 * (phi_b0_DEP_dVbs - phi_j0_DEP_dVbs) * T0 * T3 ;
732 
733          Q_sub0_dep_dVgs = - C_ESI / W_sub0 * phi_j0_DEP_dVgs ;
734          Q_sub0_dep_dVds = - C_ESI / W_sub0 * phi_j0_DEP_dVds ;
735          Q_sub0_dep_dVbs = - C_ESI / W_sub0 * (phi_j0_DEP_dVbs - Vbsc_dVbse) ;
736 
737        }
738 
739        T1 = phi_b0_DEP - phi_j0_DEP ;
740        Fn_SL_CP( T2 , T1 , 0.0 , 0.1, 2 , T0 )
741        W_b0 = sqrt(C_2ESIpq_Ndepm * (T2) ) ;
742        Fn_SU_CP( T9, W_b0, model->HSM2_tndep, 1e-8, 2 , T3 )
743        W_b0_dVgs = C_ESI / q_Ndepm / W_b0 * (phi_b0_DEP_dVgs - phi_j0_DEP_dVgs) * T0 * T3 ;
744        W_b0_dVds = C_ESI / q_Ndepm / W_b0 * (phi_b0_DEP_dVds - phi_j0_DEP_dVds) * T0 * T3 ;
745        W_b0_dVbs = C_ESI / q_Ndepm / W_b0 * (phi_b0_DEP_dVbs - phi_j0_DEP_dVbs) * T0 * T3 ;
746 
747        T1 = phi_b0_DEP - phi_s0_DEP ;
748        Fn_SL_CP( T2 , T1 , 0.0 , 0.1, 2 , T0 ) // HV_dlt=0.05
749        W_s0 = sqrt(C_2ESIpq_Ndepm * (T2) ) ;
750 
751        W_s0_dVgs = C_ESI / q_Ndepm / W_s0 * (phi_b0_DEP_dVgs - phi_s0_DEP_dVgs) * T0 ;
752        W_s0_dVds = C_ESI / q_Ndepm / W_s0 * (phi_b0_DEP_dVds - phi_s0_DEP_dVds) * T0 ;
753        W_s0_dVbs = C_ESI / q_Ndepm / W_s0 * (phi_b0_DEP_dVbs - phi_s0_DEP_dVbs) * T0 ;
754 
755        T1 = model->HSM2_tndep - T9 - W_s0 ;
756        Fn_SL_CP( W_res0 , T1 , 1.0e-20 , 1.0e-14, 2 , T0 ) // HV_dlt=1.0e-25,1.0e-18
757 
758        Qn_res0 = - W_res0 * q_Ndepm ;
759        Qn_res0_dVgs = (W_s0_dVgs + W_b0_dVgs) * q_Ndepm * T0 ;
760        Qn_res0_dVds = (W_s0_dVds + W_b0_dVds) * q_Ndepm * T0 ;
761        Qn_res0_dVbs = (W_s0_dVbs + W_b0_dVbs) * q_Ndepm * T0 ;
762 
763        if( W_bsub0 > model->HSM2_tndep && depmode !=2 ) {
764          Fn_SU_CP(T3 , phi_s0_DEP , phi_b0_DEP_lim , 0.8, 2 , T1 )
765          T3_dVgs = phi_s0_DEP_dVgs * T1 ;
766          T3_dVds = phi_s0_DEP_dVds * T1 ;
767          T3_dVbs = phi_s0_DEP_dVbs * T1 ;
768        } else {
769          Fn_SU_CP(T3 , phib_ref , phi_b0_DEP_lim , 0.8, 2 , T0 )
770          T3_dVgs = phib_ref_dVgs * T0 ;
771          T3_dVds = phib_ref_dVds * T0 ;
772          T3_dVbs = phib_ref_dVbs * T0 ;
773        }
774 
775        T4 = exp(beta * (T3 - phi_b0_DEP_lim)) ;
776        T5 = - C_QE * here->HSM2_ndepm ;
777        Qn_bac0 = T5 * T4 * T9 ;
778        Qn_bac0_dVgs = T5 * (beta * T4 * T3_dVgs * T9 + T4 * W_b0_dVgs) ;
779        Qn_bac0_dVds = T5 * (beta * T4 * T3_dVds * T9 + T4 * W_b0_dVds) ;
780        Qn_bac0_dVbs = T5 * (beta * T4 * T3_dVbs * T9 + T4 * W_b0_dVbs) ;
781 
782 
783        T1 = phi_s0_DEP - phi_b0_DEP_lim ;
784        Fn_SL_CP( T2 , T1 , 0.0, Ps_delta, 2 , T0 )
785        T2_dVgs = phi_s0_DEP_dVgs * T0 ;
786        T2_dVds = phi_s0_DEP_dVds * T0 ;
787        T2_dVbs = phi_s0_DEP_dVbs * T0 ;
788 
789        T3 = exp(beta * (T2)) ;
790        T3_dVgs = beta * T3 * T2_dVgs ;
791        T3_dVds = beta * T3 * T2_dVds ;
792        T3_dVbs = beta * T3 * T2_dVbs ;
793 
794        T4 = T3 - 1.0 - beta * T2 ;
795 
796        T4_dVgs = T3_dVgs - beta * T2_dVgs ;
797        T4_dVds = T3_dVds - beta * T2_dVds ;
798        T4_dVbs = T3_dVbs - beta * T2_dVbs ;
799 
800        T5 = sqrt(T4) ;
801        Q_n0_cur = - here->HSM2_cnst0 * T5 ;
802        Q_n0_cur_dVgs = - here->HSM2_cnst0 / 2.0 / T5 * T4_dVgs ;
803        Q_n0_cur_dVds = - here->HSM2_cnst0 / 2.0 / T5 * T4_dVds ;
804        Q_n0_cur_dVbs = - here->HSM2_cnst0 / 2.0 / T5 * T4_dVbs ;
805 
806        T4 = exp(beta * Ps_delta0) - 1.0 - beta * Ps_delta0 ;
807        T5 = sqrt(T4) ;
808        Qn_delta = here->HSM2_cnst0 * T5 ;
809 
810 
811 
812      /*-----------------------------------------------------------*
813       * Start point of phi_sL_DEP(= phi_s0_DEP + Pds) calculation.(label)
814       *-----------------*/
815 
816      /* Vdseff (begin) */
817        Vdsorg = Vds ;
818 
819        if( Vds > 1e-3 ) {
820 
821          T2 = q_Ndepm_esi / ( Cox * Cox ) ;
822          T4 = - 2.0e0 * T2 / Cox ;
823          T2_dVb = T4 * Cox_dVb ;
824          T2_dVd = T4 * Cox_dVd ;
825          T2_dVg = T4 * Cox_dVg ;
826 
827          T0 = Vgp + DEPQFN1 - beta_inv - Vbsz ;
828          T0_dVg = Vgp_dVgs ;
829          T0_dVd = Vgp_dVds - Vbsz_dVds ;
830          T0_dVb = Vgp_dVbs - Vbsz_dVbs ;
831 
832          T4 = 1.0e0 + 2.0e0 / T2 * T0 ;
833          T4_dVg = 2.0 / T2 * T0_dVg - 2.0 / T2 / T2 * T0 * T2_dVg ;
834          T4_dVd = 2.0 / T2 * T0_dVd - 2.0 / T2 / T2 * T0 * T2_dVd ;
835          T4_dVb = 2.0 / T2 * T0_dVb - 2.0 / T2 / T2 * T0 * T2_dVb ;
836 
837          Fn_SL_CP( T9 , T4 , 0 , DEPQFN_dlt, 2 , T0 )
838          T9_dVg = T4_dVg * T0 ;
839          T9_dVd = T4_dVd * T0 ;
840          T9_dVb = T4_dVb * T0 ;
841 
842          T9 +=small ;
843          T3 = sqrt( T9 ) ;
844          T3_dVg = 0.5 / T3 * T9_dVg  ;
845          T3_dVd = 0.5 / T3 * T9_dVd  ;
846          T3_dVb = 0.5 / T3 * T9_dVb  ;
847 
848          T10 = Vgp + DEPQFN1 + T2 * ( 1.0e0 - T3 ) ;
849          T10_dVb = Vgp_dVbs + T2_dVb * ( 1.0e0 - T3 ) - T2 * T3_dVb ;
850          T10_dVd = Vgp_dVds + T2_dVd * ( 1.0e0 - T3 ) - T2 * T3_dVd ;
851          T10_dVg = Vgp_dVgs + T2_dVg * ( 1.0e0 - T3 ) - T2 * T3_dVg ;
852 
853          Fn_SL_CP( T10 , T10 , DEPQFN3, 0.2, 4 , T0 )
854 //       T10 = T10 + epsm10 ;
855          T10_dVb *=T0 ;
856          T10_dVd *=T0 ;
857          T10_dVg *=T0 ;
858 
859          T1 = Vds / T10 ;
860          T2 = Fn_Pow( T1 , here->HSM2_ddlt - 1.0e0 ) ;
861          T7 = T2 * T1 ;
862          T0 = here->HSM2_ddlt * T2 / ( T10 * T10 ) ;
863          T7_dVb = T0 * ( - Vds * T10_dVb ) ;
864          T7_dVd = T0 * ( T10 - Vds * T10_dVd ) ;
865          T7_dVg = T0 * ( - Vds * T10_dVg ) ;
866 
867          T3 = 1.0 + T7 ;
868          T4 = Fn_Pow( T3 , 1.0 / here->HSM2_ddlt - 1.0 ) ;
869          T6 = T4 * T3 ;
870          T0 = T4 / here->HSM2_ddlt ;
871          T6_dVb = T0 * T7_dVb ;
872          T6_dVd = T0 * T7_dVd ;
873          T6_dVg = T0 * T7_dVg ;
874 
875          Vdseff = Vds / T6 ;
876          T0 = 1.0 / ( T6 * T6 ) ;
877          Vdseff0_dVbs = - Vds * T6_dVb * T0 ;
878          Vdseff0_dVds = ( T6 - Vds * T6_dVd ) * T0 ;
879          Vdseff0_dVgs = - Vds * T6_dVg * T0 ;
880 
881          Fn_SL_CP( Vgpp , Vgp , 0.0 , 0.5, 2 , T0 )
882          Vgpp_dVgs = T0 * Vgp_dVgs ;
883          Vgpp_dVds = T0 * Vgp_dVds ;
884          Vgpp_dVbs = T0 * Vgp_dVbs ;
885 
886          T1 = Vgpp * 0.8 ;
887          T1_dVg = Vgpp_dVgs * 0.8 ;
888          T1_dVd = Vgpp_dVds * 0.8 ;
889          T1_dVb = Vgpp_dVbs * 0.8 ;
890 
891          Fn_SU_CP3( Vds , Vdseff , Vgpp , T1, 2 , T3, T4, T5 )
892          Vdseff_dVgs = Vdseff0_dVgs * T3 + Vgpp_dVgs * T4 + T1_dVg * T5 ;
893          Vdseff_dVds = Vdseff0_dVds * T3 + Vgpp_dVds * T4 + T1_dVd * T5 ;
894          Vdseff_dVbs = Vdseff0_dVbs * T3 + Vgpp_dVbs * T4 + T1_dVb * T5 ;
895 
896        } else {
897 
898          Vdseff = Vds ;
899          Vdseff0_dVgs = 0.0 ;
900          Vdseff0_dVds = 1.0 ;
901          Vdseff0_dVbs = 0.0 ;
902 
903          Vdseff_dVgs = 0.0 ;
904          Vdseff_dVds = 1.0 ;
905          Vdseff_dVbs = 0.0 ;
906 
907        }
908        /* Vdseff (end) */
909 
910      /*---------------------------------------------------*
911       * start of phi_sL_DEP calculation. (label)
912       *--------------------------------*/
913 
914        if( Vds < 0.0e0 ) {
915 
916          phi_sL_DEP = phi_s0_DEP ;
917          phi_sL_DEP_dVgs = phi_s0_DEP_dVgs ;
918          phi_sL_DEP_dVds = phi_s0_DEP_dVds ;
919          phi_sL_DEP_dVbs = phi_s0_DEP_dVbs ;
920 
921          phi_bL_DEP = phi_b0_DEP ;
922          phi_bL_DEP_dVgs = phi_b0_DEP_dVgs ;
923          phi_bL_DEP_dVds = phi_b0_DEP_dVds ;
924          phi_bL_DEP_dVbs = phi_b0_DEP_dVbs ;
925 
926          phi_jL_DEP = phi_j0_DEP ;
927          phi_jL_DEP_dVgs = phi_j0_DEP_dVgs ;
928          phi_jL_DEP_dVds = phi_j0_DEP_dVds ;
929          phi_jL_DEP_dVbs = phi_j0_DEP_dVbs ;
930 
931          Q_subL = Q_sub0 ;
932          Q_subL_dVgs = Q_sub0_dVgs ;
933          Q_subL_dVds = Q_sub0_dVds ;
934          Q_subL_dVbs = Q_sub0_dVbs ;
935 
936          Q_nL = Q_n0 ;
937          Q_nL_dVgs = Q_n0_dVgs ;
938          Q_nL_dVds = Q_n0_dVds ;
939          Q_nL_dVbs = Q_n0_dVbs ;
940 
941 	 Q_bL_dep = Q_b0_dep ;
942 	 Q_bL_dep_dVgs = Q_b0_dep_dVgs ;
943 	 Q_bL_dep_dVds = Q_b0_dep_dVds ;
944 	 Q_bL_dep_dVbs = Q_b0_dep_dVbs ;
945 
946 	 Q_subL_dep = Q_sub0_dep ;
947 	 Q_subL_dep_dVgs = Q_sub0_dep_dVgs ;
948 	 Q_subL_dep_dVds = Q_sub0_dep_dVds ;
949 	 Q_subL_dep_dVbs = Q_sub0_dep_dVbs ;
950 
951 	 Q_sL_dep = Q_s0_dep ;
952 	 Q_sL_dep_dVgs = Q_s0_dep_dVgs ;
953 	 Q_sL_dep_dVds = Q_s0_dep_dVds ;
954 	 Q_sL_dep_dVbs = Q_s0_dep_dVbs ;
955 
956 	 Q_nL_cur = Q_n0_cur ;
957 	 Q_nL_cur_dVgs = Q_n0_cur_dVgs ;
958 	 Q_nL_cur_dVds = Q_n0_cur_dVds ;
959 	 Q_nL_cur_dVbs = Q_n0_cur_dVbs ;
960 
961        } else {
962 
963          W_bsubL = sqrt(C_2ESIpq_Ndepm * here->HSM2_nsub / (here->HSM2_nsub + here->HSM2_ndepm) * (Vds - Vbsc + Vbi_DEP)) ;
964 
965        /*---------------------------------------------------*
966         * region judgement
967         *------------------*/
968 
969         /* fully depleted case */
970          if( W_bsubL > model->HSM2_tndep ) {
971 
972            Vgp0 = Vds ;
973            W_bL = model->HSM2_tndep ;
974            phi_bL_DEP = Vds ;
975            phi_bL_DEP_lim = Vds ;
976            phi_jL_DEP = phi_bL_DEP - C_2ESIpq_Ndepm_inv * W_bL * W_bL ;
977 
978            Vgp0old = Vgp0 ;
979            phi_jL_DEP_old = phi_jL_DEP ;
980 
981            Q_bL_dep = W_bL * q_Ndepm ;
982 
983            for ( lp_s0 = 1 ; lp_s0 <= lp_se_max + 1 ; lp_s0 ++ ) {
984 
985              W_bL = sqrt(C_2ESIpq_Ndepm * (phi_bL_DEP - phi_jL_DEP) ) ;
986              Fn_SU_CP( W_bL , W_bL , model->HSM2_tndep , 1e-8, 2 , T0 )
987              W_subL = sqrt(C_2ESIpq_Nsub * (phi_jL_DEP - Vbsc + Vbi_DEP) ) ;
988 
989              Q_bL_dep = W_bL * q_Ndepm ;
990              Q_bL_dep_dPd = - C_ESI / W_bL * T0 ;
991              Q_subL_dep = - W_subL * q_Nsub ;
992              Q_subL_dep_dPd = - C_ESI / W_subL ;
993 
994              y1 = Cox * (Vgp0 - phi_bL_DEP) + Q_bL_dep + Q_subL_dep ;
995              y11 = Cox ;
996              y12 = Q_bL_dep_dPd + Q_subL_dep_dPd ;
997 
998              y2 = phi_jL_DEP - NdepmpNsub_inv1 * (NdepmpNsub * phi_bL_DEP + Vbsc - Vbi_DEP) ;
999              y21 = 0.0 ;
1000              y22 = 1.0 ;
1001 
1002              dety = y11 * y22 - y21 * y12;
1003              rev11 = (y22) / dety ;
1004              rev12 = ( - y12) / dety ;
1005              rev21 = ( - y21) / dety ;
1006              rev22 = (y11) / dety ;
1007 
1008              if( fabs( rev11 * y1 + rev12 * y2 ) > 0.5 ) {
1009                Vgp0 = Vgp0 - 0.5 * Fn_Sgn( rev11 * y1 + rev12 * y2 ) ;
1010                phi_jL_DEP = phi_jL_DEP - 0.5 * Fn_Sgn( rev21 * y1 + rev22 * y2 ) ;
1011              } else {
1012                Vgp0 = Vgp0 - ( rev11 * y1 + rev12 * y2 ) ;
1013                phi_jL_DEP = phi_jL_DEP - ( rev21 * y1 + rev22 * y2 ) ;
1014              }
1015 
1016              if( fabs(Vgp0 - Vgp0old) <= ps_conv &&
1017                  fabs(phi_jL_DEP - phi_jL_DEP_old) <= ps_conv ) lp_s0=lp_se_max + 1 ;
1018 
1019              Vgp0old = Vgp0 ;
1020              phi_jL_DEP_old = phi_jL_DEP ;
1021            }
1022            phi_jL_DEP_acc = phi_jL_DEP ;
1023 
1024            W_subL = model->HSM2_tndep * NdepmpNsub ;
1025            phi_jL_DEP = C_2ESIpq_Nsub_inv * W_subL * W_subL + Vbsc - Vbi_DEP ;
1026            phi_bL_DEP = phi_jL_DEP + C_2ESIpq_Ndepm_inv * Tn2 ;
1027            phi_sL_DEP = phi_bL_DEP ;
1028            Psbmax = phi_bL_DEP ;
1029            Vgp1 = phi_bL_DEP ;
1030            if( Vgp > Vgp0 ) {
1031              depmode = 1 ;
1032            } else if(Vgp > Vgp1 ) {
1033              depmode = 3 ;
1034            } else {
1035              depmode = 2 ;
1036            }
1037 
1038         /* else */
1039          } else {
1040            Vgp0 = Vds ;
1041            Vgp1 = Vgp0 ;
1042            Psbmax = Vgp0 ;
1043            phi_bL_DEP_lim = Vgp0 ;
1044            W_bL = W_bsubL ;
1045            W_subL = W_bL * NdepmpNsub ;
1046            phi_jL_DEP = C_2ESIpq_Nsub_inv * W_subL * W_subL + Vbsc - Vbi_DEP ;
1047            phi_bL_DEP = C_2ESIpq_Ndepm_inv * W_bL * W_bL + phi_jL_DEP ;
1048            phi_jL_DEP_acc = phi_jL_DEP ;
1049            if( Vgp > Vgp0 ) {
1050              depmode = 1 ;
1051            } else {
1052              depmode = 2 ;
1053            }
1054 
1055          }
1056 
1057          T1 = C_2ESI_q_Ndepm * ( Psbmax - ( - here->HSM2_Pb2n + Vbsc)) ;
1058          if ( T1 > 0.0 ) {
1059            vthn = - here->HSM2_Pb2n + Vbsc - sqrt(T1) / Cox ;
1060          } else {
1061            vthn = - here->HSM2_Pb2n + Vbsc ;
1062          }
1063 
1064        /*---------------------------------------------------*
1065         * initial potential phi_s0_DEP,phi_bL_DEP,phi_jL_DEP calculated.
1066         *------------------*/
1067 
1068 
1069         /* accumulation region */
1070          if( Vgp > Vgp0 ) {
1071            phi_jL_DEP = phi_jL_DEP_acc ;
1072            phi_bL_DEP = Vds ;
1073            phi_sL_DEP_ini = log(afact * Vgp * Vgp) / (beta + 2.0 / Vgp) + Vds ;
1074 
1075            if( phi_sL_DEP_ini < phi_bL_DEP_lim + ps_conv23 ) phi_sL_DEP_ini = phi_bL_DEP_lim + ps_conv23 ;
1076 
1077         /* fully depleted region */
1078          } else if( Vgp > Vgp1 ) {
1079 
1080            phi_sL_DEP_ini = phi_sL_DEP ;
1081 
1082         /* depletion & inversion */
1083 
1084          } else {
1085 
1086           /* depletion */
1087            if( Vgp > vthn ) {
1088              bfact = - 2.0 * afact * Vgp + beta ;
1089              cfact = afact * Vgp * Vgp - beta * phi_bL_DEP ;
1090              phi_bL_DEP_old = phi_bL_DEP ;
1091              phi_sL_DEP_ini = ( - bfact + sqrt(bfact * bfact - 4.0 * afact * cfact)) / 2.0 / afact ;
1092              if( phi_sL_DEP_ini > Psbmax - ps_conv23 ) phi_sL_DEP_ini = Psbmax - ps_conv23 ;
1093              W_sL = sqrt(C_2ESIpq_Ndepm * (phi_bL_DEP - phi_sL_DEP_ini) ) ;
1094              W_bL = sqrt(C_2ESIpq_Ndepm * (phi_bL_DEP - phi_jL_DEP) ) ;
1095 
1096              if( W_sL + W_bL > model->HSM2_tndep ) {
1097                for ( lp_s0 = 1 ; lp_s0 <= lp_se_max + 1 ; lp_s0 ++ ) {
1098 
1099                  y0 = W_sL + W_bL - model->HSM2_tndep ;
1100 
1101                  dydPsm = C_ESI / q_Ndepm / W_sL
1102                           + C_ESI / q_Ndepm * ( 1.0 - (here->HSM2_ndepm
1103                           / here->HSM2_nsub) / ( 1.0 + (NdepmpNsub))) / W_bL ;
1104 
1105                  if( fabs(y0 / dydPsm) > 0.5 ) {
1106                    phi_bL_DEP = phi_bL_DEP - 0.5 * Fn_Sgn(y0 / dydPsm) ;
1107                  } else {
1108                    phi_bL_DEP = phi_bL_DEP - y0 / dydPsm ;
1109                  }
1110                  if( (phi_bL_DEP - Vbsc + Vbi_DEP) < epsm10 )
1111                     phi_bL_DEP=Vbsc - Vbi_DEP + epsm10 ;
1112 
1113                  cfact = afact * Vgp * Vgp - beta * phi_bL_DEP ;
1114                  T1 = bfact * bfact - 4.0 * afact * cfact ;
1115                  if( T1 > 0.0 ) {
1116                    phi_sL_DEP_ini = ( - bfact + sqrt(T1)) / 2.0 / afact ;
1117                  } else {
1118                    phi_sL_DEP_ini = ( - bfact) / 2.0 / afact ;
1119                  }
1120 
1121                  if( phi_sL_DEP_ini > Psbmax ) phi_sL_DEP_ini = Psbmax ;
1122                  if( phi_sL_DEP_ini > phi_bL_DEP ) {
1123                    phi_sL_DEP_ini = phi_bL_DEP - ps_conv23 ;
1124                    lp_s0=lp_se_max + 1 ;
1125                  }
1126                  W_sL = sqrt(C_2ESIpq_Ndepm * (phi_bL_DEP - phi_sL_DEP_ini) ) ;
1127                  phi_jL_DEP = ( NdepmpNsub * phi_bL_DEP
1128                    + Vbsc - Vbi_DEP) / (1.0 + NdepmpNsub) ;
1129                  W_bL = sqrt(C_2ESIpq_Ndepm * (phi_bL_DEP - phi_jL_DEP) ) ;
1130 
1131                  if( fabs(phi_bL_DEP - phi_bL_DEP_old) <= 1.0e-8 ) lp_s0=lp_se_max + 1 ;
1132                  phi_bL_DEP_old = phi_bL_DEP ;
1133                }
1134              }
1135 
1136           /* inversion */
1137            } else {
1138 
1139              phi_bL_DEP = phi_b0_DEP ;
1140              phi_jL_DEP = phi_j0_DEP ;
1141              phi_sL_DEP_ini = phi_s0_DEP ;
1142 
1143            }
1144 
1145          }
1146 
1147          phi_b0_DEP_ini = phi_bL_DEP ;
1148          /*                              */
1149          /* solve poisson  at drain side */
1150          /*                              */
1151 
1152          flg_conv = 0 ;
1153 
1154         /* accumulation */
1155 
1156          phi_sL_DEP = phi_sL_DEP_ini ;
1157          phi_bL_DEP = phi_b0_DEP_ini ;
1158 
1159          phi_sL_DEP_old = phi_sL_DEP ;
1160          phi_bL_DEP_old = phi_bL_DEP ;
1161 
1162          for ( lp_s0 = 1 ; lp_s0 <= lp_se_max + 1 ; lp_s0 ++ ) {
1163 
1164            phi_jL_DEP  = NdepmpNsub_inv1 * (NdepmpNsub * phi_bL_DEP + Vbsc - Vbi_DEP) ;
1165            phi_jL_DEP_dPb = NdepmpNsub_inv1 * NdepmpNsub ;
1166 
1167            T1 = phi_bL_DEP - phi_jL_DEP ;
1168            Fn_SL_CP( T2 , T1 , 0.0 , 0.1, 2 , T7 )
1169            W_bL = sqrt(C_2ESIpq_Ndepm * (T2) ) ;
1170            Fn_SU_CP( W_bL , W_bL , model->HSM2_tndep , 1e-8, 2 , T8 )
1171            W_subL = sqrt(C_2ESIpq_Nsub * (phi_jL_DEP - Vbsc + Vbi_DEP) ) ;
1172            Q_bL_dep = W_bL * q_Ndepm ;
1173            Q_bL_dep_dPb = C_ESI / W_bL * T7 * T8 ;
1174            Q_bL_dep_dPd = - C_ESI / W_bL * T7 * T8 ;
1175            Q_subL_dep = - W_subL * q_Nsub ;
1176            Q_subL_dep_dPd = - C_ESI / W_subL ;
1177 
1178            T10 = 8.0 * q_Ndepm_esi * Tn2 ;
1179            phib_ref = (4.0 * phi_jL_DEP * phi_jL_DEP * C_ESI2 - 8.0 * phi_jL_DEP * C_ESI2 * phi_sL_DEP
1180                     + 4.0 * C_ESI2 * phi_sL_DEP * phi_sL_DEP
1181                     + 4.0 * phi_jL_DEP * q_Ndepm_esi * Tn2
1182                     + 4.0 * phi_sL_DEP * q_Ndepm_esi * Tn2
1183                     + Ndepm2 * C_QE2 * Tn2 * Tn2  ) / T10 ;
1184            phib_ref_dPs = ( - 8.0 * phi_jL_DEP * C_ESI2 + 8.0 * C_ESI2 * phi_sL_DEP
1185                         + 4.0 * q_Ndepm_esi * Tn2 ) / T10 ;
1186            phib_ref_dPd = (   8.0 * phi_jL_DEP * C_ESI2 - 8.0 * C_ESI2 * phi_sL_DEP
1187                         + 4.0 * q_Ndepm_esi * Tn2 ) / T10 ;
1188 
1189            T1 = beta * (phi_sL_DEP - phi_bL_DEP) ;
1190            T2 = exp(T1) ;
1191            if( phi_sL_DEP >= phi_bL_DEP ) {
1192              Q_sL = - here->HSM2_cnst0 * sqrt(T2 - 1.0 - T1 + 1e-15) ;
1193              Q_sL_dPs = 0.5 * here->HSM2_cnst0 * here->HSM2_cnst0 / Q_sL * (beta * T2 - beta) ;
1194              Q_sL_dPb = - Q_sL_dPs ;
1195            } else {
1196              T3 = exp( - beta * (phi_sL_DEP - Vbsc)) ;
1197              T4 = exp( - beta * (phi_bL_DEP - Vbsc)) ;
1198              Q_sL = here->HSM2_cnst0 * sqrt(T2 - 1.0 - T1 + here->HSM2_cnst1 * (T3 - T4) + 1e-15) ;
1199              T5 = 0.5 * here->HSM2_cnst0 * here->HSM2_cnst0 / Q_sL ;
1200              Q_sL_dPs = T5 * (beta * T2 - beta + here->HSM2_cnst1 * ( - beta * T3) ) ;
1201              Q_sL_dPb = T5 * ( - beta * T2 + beta + here->HSM2_cnst1 * beta * T4 ) ;
1202            }
1203 
1204            Fn_SU_CP2( T1 , phib_ref , phi_bL_DEP_lim , sm_delta, 4 , T9, T11 )
1205 
1206            y1 = phi_bL_DEP - T1 ;
1207            y11 = - phib_ref_dPs * T9 ;
1208            y12 = 1.0 - phib_ref_dPd * phi_jL_DEP_dPb * T9 ;
1209 
1210            y2 = Cox * (Vgp - phi_sL_DEP) + Q_sL + Q_bL_dep + Q_subL_dep ;
1211            y21 = - Cox + Q_sL_dPs ;
1212            y22 = Q_sL_dPb + Q_bL_dep_dPb + Q_bL_dep_dPd * phi_jL_DEP_dPb + Q_subL_dep_dPd * phi_jL_DEP_dPb ;
1213 
1214            dety = y11 * y22 - y21 * y12;
1215            rev11 = (y22) / dety ;
1216            rev12 = ( - y12) / dety ;
1217            rev21 = ( - y21) / dety ;
1218            rev22 = (y11) / dety ;
1219            if( fabs( rev21 * y1 + rev22 * y2 ) > 0.2 ) { // HV_dlt=0.5
1220              phi_sL_DEP = phi_sL_DEP - 0.2 * Fn_Sgn( rev11 * y1 + rev12 * y2 ) ;
1221              phi_bL_DEP = phi_bL_DEP - 0.2 * Fn_Sgn( rev21 * y1 + rev22 * y2 ) ;
1222            } else {
1223              phi_sL_DEP = phi_sL_DEP - ( rev11 * y1 + rev12 * y2 ) ;
1224              phi_bL_DEP = phi_bL_DEP - ( rev21 * y1 + rev22 * y2 ) ;
1225            }
1226 
1227            if( fabs(phi_sL_DEP - phi_sL_DEP_old) <= ps_conv &&  fabs(phi_bL_DEP - phi_bL_DEP_old) <= ps_conv ) {
1228              lp_s0=lp_se_max + 1 ;
1229              flg_conv = 1 ;
1230            }
1231 
1232            phi_sL_DEP_old = phi_sL_DEP ;
1233            phi_bL_DEP_old = phi_bL_DEP ;
1234 
1235          }
1236          if( flg_conv == 0 ) {
1237            printf( "*** warning(HiSIM(%s)): Went Over Iteration Maximum(Psl)\n",model->HSM2modName ) ;
1238            printf( " Vbse   = %7.3f Vdse = %7.3f Vgse = %7.3f\n" ,Vbse , Vdse , Vgse ) ;
1239          }
1240 
1241         /* caluculate derivative */
1242 
1243          y1_dVgs = - Vdseff_dVgs * T11 ;
1244          y1_dVds = - Vdseff_dVds * T11 ;
1245          y1_dVbs = - (8.0 * phi_jL_DEP * C_ESI2 - 8.0 * C_ESI2 * phi_sL_DEP
1246                   + 4.0 * q_Ndepm_esi * Tn2) / T10
1247                   * T9 * NdepmpNsub_inv1 * Vbsc_dVbse - Vdseff_dVbs * T11 ;
1248 
1249          Q_bL_dep_dVbs = - C_ESI / W_bL * T7 * T8 * NdepmpNsub_inv1 * Vbsc_dVbse ;
1250 
1251          Q_subL_dep_dVbs = - C_ESI / W_subL * (NdepmpNsub_inv1 * Vbsc_dVbse - Vbsc_dVbse) ;
1252 
1253          T1 = beta * (phi_sL_DEP - phi_bL_DEP) ;
1254          T2 = exp(T1) ;
1255          if( phi_sL_DEP >= phi_bL_DEP ) {
1256            Q_sL_dVbs = 0.0 ;
1257          } else {
1258            T3 = exp( - beta * (phi_sL_DEP - Vbsc)) ;
1259            T4 = exp( - beta * (phi_bL_DEP - Vbsc)) ;
1260            T5 = sqrt(T2 - 1.0 - T1 + 1e-15 + here->HSM2_cnst1 * (T3 - T4)) ;
1261            Q_sL_dVbs = here->HSM2_cnst0 / 2.0 / T5 *
1262                   (here->HSM2_cnst1 * (beta * T3 * Vbsc_dVbse - beta * T4 * Vbsc_dVbse) ) ;
1263          }
1264 
1265          y2_dVgs = Cox_dVg * (Vgp - phi_sL_DEP) + Cox * Vgp_dVgs ;
1266          y2_dVds = Cox_dVd * (Vgp - phi_sL_DEP) + Cox * Vgp_dVds ;
1267          y2_dVbs = Cox_dVb * (Vgp - phi_sL_DEP) + Cox * Vgp_dVbs + Q_sL_dVbs + Q_bL_dep_dVbs + Q_subL_dep_dVbs ;
1268 
1269          phi_sL_DEP_dVgs = - ( rev11 * y1_dVgs + rev12 * y2_dVgs ) ;
1270          phi_sL_DEP_dVds = - ( rev11 * y1_dVds + rev12 * y2_dVds ) ;
1271          phi_sL_DEP_dVbs = - ( rev11 * y1_dVbs + rev12 * y2_dVbs ) ;
1272 
1273          phi_bL_DEP_dVgs = - ( rev21 * y1_dVgs + rev22 * y2_dVgs ) ;
1274          phi_bL_DEP_dVds = - ( rev21 * y1_dVds + rev22 * y2_dVds ) ;
1275          phi_bL_DEP_dVbs = - ( rev21 * y1_dVbs + rev22 * y2_dVbs ) ;
1276 
1277          if( W_bsubL > model->HSM2_tndep && depmode !=2 ) {
1278            Fn_SU_CP2(phi_bL_DEP , phi_bL_DEP , phi_sL_DEP , 0.04, 2 , T1, T2 ) // HV_dlt=0.02
1279            phi_bL_DEP_dVgs = phi_bL_DEP_dVgs * T1 + phi_sL_DEP_dVgs * T2 ;
1280            phi_bL_DEP_dVds = phi_bL_DEP_dVds * T1 + phi_sL_DEP_dVds * T2 ;
1281            phi_bL_DEP_dVbs = phi_bL_DEP_dVbs * T1 + phi_sL_DEP_dVbs * T2 ;
1282          }
1283 
1284          phi_jL_DEP  = NdepmpNsub_inv1 * (NdepmpNsub * phi_bL_DEP + Vbsc - Vbi_DEP) ;
1285          phi_jL_DEP_dVgs = NdepmpNsub_inv1 * NdepmpNsub * phi_bL_DEP_dVgs ;
1286          phi_jL_DEP_dVds = NdepmpNsub_inv1 * NdepmpNsub * phi_bL_DEP_dVds ;
1287          phi_jL_DEP_dVbs = NdepmpNsub_inv1 * NdepmpNsub * phi_bL_DEP_dVbs + NdepmpNsub_inv1 * Vbsc_dVbse  ;
1288 
1289          phib_ref = (4.0 * phi_jL_DEP * phi_jL_DEP * C_ESI2 - 8.0 * phi_jL_DEP * C_ESI2 * phi_sL_DEP
1290                   + 4.0 * C_ESI2 * phi_sL_DEP * phi_sL_DEP
1291                   + 4.0 * phi_jL_DEP * q_Ndepm_esi * Tn2
1292                   + 4.0 * phi_sL_DEP * q_Ndepm_esi * Tn2
1293                   + Ndepm2 * C_QE2 * Tn2 * Tn2 ) / T10 ;
1294 
1295          phib_ref_dVgs = ( 8.0 * phi_jL_DEP * phi_jL_DEP_dVgs * C_ESI2 - 8.0 * phi_jL_DEP_dVgs * C_ESI2 * phi_sL_DEP
1296                          - 8.0 * phi_jL_DEP * C_ESI2 * phi_sL_DEP_dVgs + 8.0 * C_ESI2 * phi_sL_DEP * phi_sL_DEP_dVgs
1297                        + 4.0 * phi_jL_DEP_dVgs * q_Ndepm_esi * Tn2
1298                        + 4.0 * phi_sL_DEP_dVgs * q_Ndepm_esi * Tn2 ) / T10 ;
1299          phib_ref_dVds = ( 8.0 * phi_jL_DEP * phi_jL_DEP_dVds * C_ESI2 - 8.0 * phi_jL_DEP_dVds * C_ESI2 * phi_sL_DEP
1300                          - 8.0 * phi_jL_DEP * C_ESI2 * phi_sL_DEP_dVds + 8.0 * C_ESI2 * phi_sL_DEP * phi_sL_DEP_dVds
1301                        + 4.0 * phi_jL_DEP_dVds * q_Ndepm_esi * Tn2
1302                        + 4.0 * phi_sL_DEP_dVds * q_Ndepm_esi * Tn2 ) / T10 ;
1303          phib_ref_dVbs = ( 8.0 * phi_jL_DEP * phi_jL_DEP_dVbs * C_ESI2 - 8.0 * phi_jL_DEP_dVbs * C_ESI2 * phi_sL_DEP
1304                          - 8.0 * phi_jL_DEP * C_ESI2 * phi_sL_DEP_dVbs + 8.0 * C_ESI2 * phi_sL_DEP * phi_sL_DEP_dVbs
1305                        + 4.0 * phi_jL_DEP_dVbs * q_Ndepm_esi * Tn2
1306                        + 4.0 * phi_sL_DEP_dVbs * q_Ndepm_esi * Tn2 ) / T10 ;
1307 
1308          T1 = beta * (phi_sL_DEP - phi_bL_DEP) ;
1309          T1_dVgs = beta * (phi_sL_DEP_dVgs - phi_bL_DEP_dVgs) ;
1310          T1_dVds = beta * (phi_sL_DEP_dVds - phi_bL_DEP_dVds) ;
1311          T1_dVbs = beta * (phi_sL_DEP_dVbs - phi_bL_DEP_dVbs) ;
1312 
1313          T2 = exp(T1) ;
1314          T2_dVgs = T1_dVgs * T2 ;
1315          T2_dVds = T1_dVds * T2 ;
1316          T2_dVbs = T1_dVbs * T2 ;
1317 
1318          if( phi_sL_DEP >= phi_bL_DEP ) {
1319            T3 = sqrt(T2 - 1.0e0 - T1 + 1e-15 ) ;
1320            T3_dVgs = (T2_dVgs - T1_dVgs) / 2.0 / T3 ;
1321            T3_dVds = (T2_dVds - T1_dVds) / 2.0 / T3 ;
1322            T3_dVbs = (T2_dVbs - T1_dVbs) / 2.0 / T3 ;
1323 
1324            Q_sL = - here->HSM2_cnst0 * T3 ;
1325 
1326            Q_sL_dep = 0.0 ;
1327            Q_subL = 0.0 ;
1328 
1329            W_bL = sqrt(C_2ESIpq_Ndepm * (phi_bL_DEP - phi_jL_DEP) ) ;
1330            Fn_SU_CP( T9 , W_bL , model->HSM2_tndep, 1e-8, 2 , T4 )
1331 
1332            W_subL = sqrt(C_2ESIpq_Nsub * (phi_jL_DEP - Vbsc + Vbi_DEP) ) ;
1333            Q_bL_dep = T9 * q_Ndepm ;
1334            Q_subL_dep = - W_subL * q_Nsub ;
1335 
1336           /* derivative */
1337            Q_sL_dVgs = - here->HSM2_cnst0 * T3_dVgs ;
1338            Q_sL_dVds = - here->HSM2_cnst0 * T3_dVds ;
1339            Q_sL_dVbs = - here->HSM2_cnst0 * T3_dVbs ;
1340 
1341            Q_nL = Q_sL ;
1342            Q_nL_dVgs = Q_sL_dVgs ;
1343            Q_nL_dVds = Q_sL_dVds ;
1344            Q_nL_dVbs = Q_sL_dVbs ;
1345 
1346            Q_bL_dep_dVgs = C_ESI / W_bL * (phi_bL_DEP_dVgs - phi_jL_DEP_dVgs) * T4 ;
1347            Q_bL_dep_dVds = C_ESI / W_bL * (phi_bL_DEP_dVds - phi_jL_DEP_dVds) * T4 ;
1348            Q_bL_dep_dVbs = C_ESI / W_bL * (phi_bL_DEP_dVbs - phi_jL_DEP_dVbs) * T4 ;
1349 
1350            Q_subL_dep_dVgs = - C_ESI / W_subL * phi_jL_DEP_dVgs ;
1351            Q_subL_dep_dVds = - C_ESI / W_subL * phi_jL_DEP_dVds ;
1352            Q_subL_dep_dVbs = - C_ESI / W_subL * (phi_jL_DEP_dVbs - Vbsc_dVbse) ;
1353 
1354            Q_subL_dVgs = 0.0 ;
1355            Q_subL_dVds = 0.0 ;
1356            Q_subL_dVbs = 0.0 ;
1357 
1358            Q_sL_dep_dVgs = 0.0 ;
1359            Q_sL_dep_dVds = 0.0 ;
1360            Q_sL_dep_dVbs = 0.0 ;
1361 
1362          } else {
1363 
1364            T3 = exp( - beta * (phi_sL_DEP - Vbsc)) ;
1365            T4 = exp( - beta * (phi_bL_DEP - Vbsc)) ;
1366            T5 = sqrt(T2 - 1.0 - T1 + here->HSM2_cnst1 * (T3 - T4) + 1e-15) ;
1367            Q_sL = here->HSM2_cnst0 * T5 ;
1368 
1369            T3_dVgs = - beta * T3 * phi_sL_DEP_dVgs ;
1370            T3_dVds = - beta * T3 * phi_sL_DEP_dVds ;
1371            T3_dVbs = - beta * T3 * (phi_sL_DEP_dVbs - Vbsc_dVbse) ;
1372 
1373            T4_dVgs = - beta * T4 * phi_bL_DEP_dVgs ;
1374            T4_dVds = - beta * T4 * phi_bL_DEP_dVds ;
1375            T4_dVbs = - beta * T4 * (phi_bL_DEP_dVbs - Vbsc_dVbse) ;
1376 
1377            T5_dVgs = (T2_dVgs - T1_dVgs + here->HSM2_cnst1 * (T3_dVgs - T4_dVgs)) / 2.0 / T5 ;
1378            T5_dVds = (T2_dVds - T1_dVds + here->HSM2_cnst1 * (T3_dVds - T4_dVds)) / 2.0 / T5 ;
1379            T5_dVbs = (T2_dVbs - T1_dVbs + here->HSM2_cnst1 * (T3_dVbs - T4_dVbs)) / 2.0 / T5 ;
1380 
1381            Q_sL_dVgs = here->HSM2_cnst0 * T5_dVgs ;
1382            Q_sL_dVds = here->HSM2_cnst0 * T5_dVds ;
1383            Q_sL_dVbs = here->HSM2_cnst0 * T5_dVbs ;
1384 
1385            if( W_bsubL > model->HSM2_tndep && depmode !=2 ) {
1386              Q_subL = 0.0 ;
1387              Q_sL_dep = 0.0 ;
1388 
1389              Q_subL_dVgs = 0.0 ;
1390              Q_subL_dVds = 0.0 ;
1391              Q_subL_dVbs = 0.0 ;
1392 
1393              Q_sL_dep_dVgs = 0.0 ;
1394              Q_sL_dep_dVds = 0.0 ;
1395              Q_sL_dep_dVbs = 0.0 ;
1396 
1397            } else {
1398              T3 = exp( - beta * (phi_sL_DEP - Vbsc)) ;
1399              T4 = exp( - beta * (phi_bL_DEP - Vbsc)) ;
1400              T5 = sqrt( - T1 + here->HSM2_cnst1 * (T3 - T4)) ;
1401              Q_subL = here->HSM2_cnst0 * T5 - here->HSM2_cnst0 * sqrt( - T1)  ;
1402              T6 = sqrt(T2 - 1.0e0 - T1 + 1e-15) ;
1403              Q_sL_dep = here->HSM2_cnst0 * T6 ;
1404 
1405              Q_subL_dVgs = here->HSM2_cnst0 / 2.0 / T5 * ( - beta * (phi_sL_DEP_dVgs - phi_bL_DEP_dVgs)
1406                 + here->HSM2_cnst1 * ( - beta * T3 * phi_sL_DEP_dVgs + beta * T4 * phi_bL_DEP_dVgs))
1407                 - here->HSM2_cnst0 / 2.0 / sqrt( - T1) * ( - beta * (phi_sL_DEP_dVgs - phi_bL_DEP_dVgs)) ;
1408              Q_subL_dVds = here->HSM2_cnst0 / 2.0 / T5 * ( - beta * (phi_sL_DEP_dVds - phi_bL_DEP_dVds)
1409                 + here->HSM2_cnst1 * ( - beta * T3 * phi_sL_DEP_dVds + beta * T4 * phi_bL_DEP_dVds))
1410                 - here->HSM2_cnst0 / 2.0 / sqrt( - T1) * ( - beta * (phi_sL_DEP_dVds - phi_bL_DEP_dVds)) ;
1411              Q_subL_dVbs = here->HSM2_cnst0 / 2.0 / T5 * ( - beta * (phi_sL_DEP_dVbs - phi_bL_DEP_dVbs)
1412                 + here->HSM2_cnst1 * ( - beta * T3 * (phi_sL_DEP_dVbs - Vbsc_dVbse) + beta * T4 * (phi_bL_DEP_dVbs - Vbsc_dVbse)))
1413                 - here->HSM2_cnst0 / 2.0 / sqrt( - T1) * ( - beta * (phi_sL_DEP_dVbs - phi_bL_DEP_dVbs)) ;
1414 
1415              Q_sL_dep_dVgs = here->HSM2_cnst0 / 2.0 / T6 * beta * (phi_sL_DEP_dVgs - phi_bL_DEP_dVgs) * (T2 - 1) ;
1416              Q_sL_dep_dVds = here->HSM2_cnst0 / 2.0 / T6 * beta * (phi_sL_DEP_dVds - phi_bL_DEP_dVds) * (T2 - 1) ;
1417              Q_sL_dep_dVbs = here->HSM2_cnst0 / 2.0 / T6 * beta * (phi_sL_DEP_dVbs - phi_bL_DEP_dVbs) * (T2 - 1) ;
1418 
1419            }
1420 
1421            Q_nL = 0.0 ;
1422            Q_nL_dVgs = 0.0 ;
1423            Q_nL_dVds = 0.0 ;
1424            Q_nL_dVbs = 0.0 ;
1425 
1426 ///        Qg = Cox * (Vgp - phi_sL_DEP) ;
1427 
1428            T1 = phi_bL_DEP - phi_jL_DEP ;
1429            Fn_SL_CP( T2 , T1 , 0.0 , 0.1, 2 , T0 )
1430            W_bL = sqrt(C_2ESIpq_Ndepm * (T2) ) ;
1431            Fn_SU_CP( T9 , W_bL , model->HSM2_tndep, 1e-8, 2 , T3 )
1432            W_subL = sqrt(C_2ESIpq_Nsub * (phi_jL_DEP - Vbsc + Vbi_DEP) ) ;
1433            Q_bL_dep = T9 * q_Ndepm ;
1434            Q_subL_dep = - W_subL * q_Nsub ;
1435 
1436            Q_bL_dep_dVgs = C_ESI / W_bL * (phi_bL_DEP_dVgs - phi_jL_DEP_dVgs) * T0 * T3 ;
1437            Q_bL_dep_dVds = C_ESI / W_bL * (phi_bL_DEP_dVds - phi_jL_DEP_dVds) * T0 * T3 ;
1438            Q_bL_dep_dVbs = C_ESI / W_bL * (phi_bL_DEP_dVbs - phi_jL_DEP_dVbs) * T0 * T3 ;
1439 
1440            Q_subL_dep_dVgs = - C_ESI / W_subL * phi_jL_DEP_dVgs ;
1441            Q_subL_dep_dVds = - C_ESI / W_subL * phi_jL_DEP_dVds ;
1442            Q_subL_dep_dVbs = - C_ESI / W_subL * (phi_jL_DEP_dVbs - Vbsc_dVbse) ;
1443 
1444          }
1445 
1446 
1447          T1 = phi_sL_DEP - phi_bL_DEP_lim ;
1448          Fn_SL_CP( T2 , T1 , 0.0, Ps_delta, 2 , T0 )
1449          T2_dVgs = (phi_sL_DEP_dVgs - Vdseff_dVgs) * T0 ;
1450          T2_dVds = (phi_sL_DEP_dVds - Vdseff_dVds) * T0 ;
1451          T2_dVbs = (phi_sL_DEP_dVbs - Vdseff_dVbs) * T0 ;
1452 
1453          T3 = exp(beta * (T2)) ;
1454          T3_dVgs = beta * T3 * T2_dVgs ;
1455          T3_dVds = beta * T3 * T2_dVds ;
1456          T3_dVbs = beta * T3 * T2_dVbs ;
1457 
1458          T4 = T3 - 1.0 - beta * T2 ;
1459 
1460          T4_dVgs = T3_dVgs - beta * T2_dVgs ;
1461          T4_dVds = T3_dVds - beta * T2_dVds ;
1462          T4_dVbs = T3_dVbs - beta * T2_dVbs ;
1463 
1464          T5 = sqrt(T4) ;
1465          Q_nL_cur = - here->HSM2_cnst0 * T5 ;
1466          Q_nL_cur_dVgs = - here->HSM2_cnst0 / 2.0 / T5 * T4_dVgs ;
1467          Q_nL_cur_dVds = - here->HSM2_cnst0 / 2.0 / T5 * T4_dVds ;
1468          Q_nL_cur_dVbs = - here->HSM2_cnst0 / 2.0 / T5 * T4_dVbs ;
1469 
1470        }
1471 
1472      /*---------------------------------------------------*
1473       * Assign Pds.
1474       *-----------------*/
1475        Ps0 = phi_s0_DEP ;
1476        Psl = phi_sL_DEP ;
1477        Pds = phi_sL_DEP - phi_s0_DEP + epsm10 ;
1478 
1479        Pds_dVgs = phi_sL_DEP_dVgs - phi_s0_DEP_dVgs ;
1480        Pds_dVds = phi_sL_DEP_dVds - phi_s0_DEP_dVds ;
1481        Pds_dVbs = phi_sL_DEP_dVbs - phi_s0_DEP_dVbs ;
1482 
1483        if( Pds < 0.0 ) { // take care of numerical noise //
1484          Pds = 0.0 ;
1485          Pds_dVgs = 0.0 ;
1486          Pds_dVds = 0.0 ;
1487          Pds_dVbs = 0.0 ;
1488 
1489          phi_sL_DEP = phi_s0_DEP ;
1490          phi_sL_DEP_dVgs = phi_s0_DEP_dVgs ;
1491          phi_sL_DEP_dVds = phi_s0_DEP_dVds ;
1492          phi_sL_DEP_dVbs = phi_s0_DEP_dVbs ;
1493 
1494          Idd = 0.0 ;
1495          Idd_dVgs = 0.0 ;
1496          Idd_dVds = 0.0 ;
1497          Idd_dVbs = 0.0 ;
1498 
1499        } else {
1500 
1501          T1 = - (Q_s0 + Q_sL) ;
1502          T1_dVgs = - (Q_s0_dVgs + Q_sL_dVgs) ;
1503          T1_dVds = - (Q_s0_dVds + Q_sL_dVds) ;
1504          T1_dVbs = - (Q_s0_dVbs + Q_sL_dVbs) ;
1505 
1506          Fn_SL_CP3( Qn_drift , T1 , 0.0, Qn_delta , 2 , T3, T4, T5 )
1507          Qn_drift_dVgs = T1_dVgs * T3 ;
1508          Qn_drift_dVds = T1_dVds * T3 ;
1509          Qn_drift_dVbs = T1_dVbs * T3 ;
1510 
1511          Idd_drift =  beta * Qn_drift / 2.0 * Pds ;
1512          Idd_diffu = - ( - Q_nL_cur + Q_n0_cur);
1513 
1514          Idd = Idd_drift + Idd_diffu ;
1515          Idd_dVgs = beta * Qn_drift_dVgs / 2.0 * Pds + beta * Qn_drift / 2.0 * Pds_dVgs
1516                     - ( - Q_nL_cur_dVgs + Q_n0_cur_dVgs ) ;
1517          Idd_dVds = beta * Qn_drift_dVds / 2.0 * Pds + beta * Qn_drift / 2.0 * Pds_dVds
1518                     - ( - Q_nL_cur_dVds + Q_n0_cur_dVds ) ;
1519          Idd_dVbs = beta * Qn_drift_dVbs / 2.0 * Pds + beta * Qn_drift / 2.0 * Pds_dVbs
1520                     - ( - Q_nL_cur_dVbs + Q_n0_cur_dVbs ) ;
1521 
1522        }
1523 
1524 
1525        Qiu = - Q_n0_cur ;
1526        Qiu_dVgs = - Q_n0_cur_dVgs ;
1527        Qiu_dVds = - Q_n0_cur_dVds ;
1528        Qiu_dVbs = - Q_n0_cur_dVbs ;
1529 
1530        Lch = Leff ;
1531 
1532       /*-----------------------------------------------------------*
1533        * Muun : universal mobility.  (CGS unit)
1534        *-----------------*/
1535 
1536        T2 = here->HSM2_ninv_o_esi / C_m2cm  ;
1537 
1538        T0 = here->HSM2_ninvd ;
1539        T3 = sqrt(Pds*Pds + model->HSM2_vzadd0) ;
1540        Pdsz = T3 - sqrt(model->HSM2_vzadd0) ;
1541        Pdsz_dVbs = Pds_dVbs * Pds / T3  ;
1542        Pdsz_dVds = Pds_dVds * Pds / T3  ;
1543        Pdsz_dVgs = Pds_dVgs * Pds / T3  ;
1544 
1545        T4 = 1.0 + ( Pdsz ) * T0 ;
1546        T4_dVb = ( Pdsz_dVbs ) * T0 ;
1547        T4_dVd = ( Pdsz_dVds ) * T0 ;
1548        T4_dVg = ( Pdsz_dVgs ) * T0 ;
1549 
1550        T5     = T2 * Qiu ;
1551        T5_dVb = T2 * Qiu_dVbs ;
1552        T5_dVd = T2 * Qiu_dVds ;
1553        T5_dVg = T2 * Qiu_dVgs ;
1554 
1555        T3     = T5 / T4 ;
1556        T3_dVb = ( - T4_dVb * T5 + T4 * T5_dVb ) / T4 / T4 ;
1557        T3_dVd = ( - T4_dVd * T5 + T4 * T5_dVd ) / T4 / T4 ;
1558        T3_dVg = ( - T4_dVg * T5 + T4 * T5_dVg ) / T4 / T4 ;
1559 
1560        Eeff = T3 ;
1561        Eeff_dVbs = T3_dVb ;
1562        Eeff_dVds = T3_dVd ;
1563        Eeff_dVgs = T3_dVg ;
1564 
1565        T5 = Fn_Pow( Eeff , model->HSM2_mueph0 - 1.0e0 ) ;
1566        T8 = T5 * Eeff ;
1567        T7 = Fn_Pow( Eeff , here->HSM2_muesr - 1.0e0 ) ;
1568        T6 = T7 * Eeff ;
1569 
1570        T9 = C_QE * C_m2cm_p2 ;
1571        Rns = Qiu / T9 ;
1572 
1573        T1 = 1.0e0 / ( here->HSM2_muecb0 + here->HSM2_muecb1 * Rns / 1.0e11 )
1574          + here->HSM2_mphn0 * T8 + T6 / pParam->HSM2_muesr1 ;
1575 
1576        Muun = 1.0e0 / T1 ;
1577 
1578        T1 = 1.0e0 / ( T1 * T1 ) ;
1579        T2 = here->HSM2_muecb0 + here->HSM2_muecb1 * Rns / 1.0e11 ;
1580        T2 = 1.0e0 / ( T2 * T2 ) ;
1581        /*  here->HSM2_mphn1 = here->HSM2_mphn0 * model->HSM2_mueph0 */
1582        T3 = here->HSM2_mphn1 * T5 ;
1583        T4 = here->HSM2_muesr * T7 / pParam->HSM2_muesr1 ;
1584        T5 = - 1.0e-11 * here->HSM2_muecb1 / C_QE * T2 / C_m2cm_p2 ;
1585        Muun_dVbs = - ( T5 * Qiu_dVbs
1586                     + Eeff_dVbs * T3 + Eeff_dVbs * T4 ) * T1 ;
1587        Muun_dVds = - ( T5 * Qiu_dVds
1588                     + Eeff_dVds * T3 + Eeff_dVds * T4 ) * T1 ;
1589        Muun_dVgs = - ( T5 * Qiu_dVgs
1590                     + Eeff_dVgs * T3 + Eeff_dVgs * T4 ) * T1 ;
1591 
1592       /*  Change to MKS unit */
1593        Muun /=C_m2cm_p2 ;
1594        Muun_dVbs /=C_m2cm_p2 ;
1595        Muun_dVds /=C_m2cm_p2 ;
1596        Muun_dVgs /=C_m2cm_p2 ;
1597 
1598       /*-----------------------------------------------------------*
1599        * Mu : mobility
1600        *-----------------*/
1601        T2 = beta * (Qiu + small) * Lch ;
1602 
1603        T1 = 1.0e0 / T2 ;
1604        T3 = T1 * T1 ;
1605        T4 = - beta * T3 ;
1606        T5 = T4 * Lch ;
1607        T6 = T4 * (Qiu + small) ;
1608        T1_dVb = ( T5 * Qiu_dVbs ) ;
1609        T1_dVd = ( T5 * Qiu_dVds ) ;
1610        T1_dVg = ( T5 * Qiu_dVgs ) ;
1611 
1612        TY = Idd * T1 ;
1613        TY_dVbs = Idd_dVbs * T1 + Idd * T1_dVb ;
1614        TY_dVds = Idd_dVds * T1 + Idd * T1_dVd ;
1615        TY_dVgs = Idd_dVgs * T1 + Idd * T1_dVg ;
1616 
1617        T2 = 0.2 * Vmax / Muun ;
1618        T3 = - T2 / Muun ;
1619        T2_dVb = T3 * Muun_dVbs ;
1620        T2_dVd = T3 * Muun_dVds ;
1621        T2_dVg = T3 * Muun_dVgs ;
1622 
1623        Ey = sqrt( TY * TY + T2 * T2 ) ;
1624        T4 = 1.0 / Ey ;
1625        Ey_dVbs = T4 * ( TY * TY_dVbs + T2 * T2_dVb ) ;
1626        Ey_dVds = T4 * ( TY * TY_dVds + T2 * T2_dVd ) ;
1627        Ey_dVgs = T4 * ( TY * TY_dVgs + T2 * T2_dVg ) ;
1628 
1629        Em = Muun * Ey ;
1630        Em_dVbs = Muun_dVbs * Ey + Muun * Ey_dVbs ;
1631        Em_dVds = Muun_dVds * Ey + Muun * Ey_dVds ;
1632        Em_dVgs = Muun_dVgs * Ey + Muun * Ey_dVgs ;
1633 
1634        T1  = Em / Vmax ;
1635 
1636        Ey_suf = Ey ;
1637        Ey_suf_dVgs = Ey_dVgs ;
1638        Ey_suf_dVds = Ey_dVds ;
1639        Ey_suf_dVbs = Ey_dVbs ;
1640 
1641       /* note: model->HSM2_bb = 2 (electron) ;1 (hole) */
1642        if ( 1.0e0 - epsm10 <= model->HSM2_bb && model->HSM2_bb <= 1.0e0 + epsm10 ) {
1643          T3 = 1.0e0 ;
1644        } else if ( 2.0e0 - epsm10 <= model->HSM2_bb && model->HSM2_bb <= 2.0e0 + epsm10 ) {
1645          T3 = T1 ;
1646        } else {
1647          T3 = Fn_Pow( T1 , model->HSM2_bb - 1.0e0 ) ;
1648        }
1649        T2 = T1 * T3 ;
1650        T4 = 1.0e0 + T2 ;
1651 
1652        if ( 1.0e0 - epsm10 <= model->HSM2_bb && model->HSM2_bb <= 1.0e0 + epsm10 ) {
1653          T5 = 1.0 / T4 ;
1654          T6 = T5 / T4 ;
1655        } else if ( 2.0e0 - epsm10 <= model->HSM2_bb && model->HSM2_bb <= 2.0e0 + epsm10 ) {
1656          T5 = 1.0 / sqrt( T4 ) ;
1657          T6 = T5 / T4 ;
1658        } else {
1659          T6 = Fn_Pow( T4 , ( - 1.0e0 / model->HSM2_bb - 1.0e0 ) ) ;
1660          T5 = T4 * T6 ;
1661        }
1662 
1663        T7 = Muun / Vmax * T6 * T3 ;
1664 
1665        Mu = Muun * T5 ;
1666 
1667        Mu_dVbs = Muun_dVbs * T5 - T7 * Em_dVbs ;
1668        Mu_dVds = Muun_dVds * T5 - T7 * Em_dVds ;
1669        Mu_dVgs = Muun_dVgs * T5 - T7 * Em_dVgs ;
1670 
1671      //-----------------------------------------------------------*
1672      //*  resistor region current.  (CGS unit)
1673      //*-----------------//
1674 
1675        if( Vdsorg > 1e-3 ) {
1676 
1677          T2 = q_Ndepm_esi / ( Cox * Cox ) ;
1678          T4 = - 2.0e0 * T2 / Cox ;
1679          T2_dVb = T4 * Cox_dVb ;
1680          T2_dVd = T4 * Cox_dVd ;
1681          T2_dVg = T4 * Cox_dVg ;
1682 
1683          T0 = Vgp + here->HSM2_depvdsef1 - beta_inv - Vbsz ;
1684          T0_dVg = Vgp_dVgs ;
1685          T0_dVd = Vgp_dVds - Vbsz_dVds ;
1686          T0_dVb = Vgp_dVbs - Vbsz_dVbs ;
1687 
1688          T4 = 1.0e0 + 2.0e0 / T2 * T0 ;
1689          T4_dVg = 2.0 / T2 * T0_dVg - 2.0 / T2 / T2 * T0 * T2_dVg ;
1690          T4_dVd = 2.0 / T2 * T0_dVd - 2.0 / T2 / T2 * T0 * T2_dVd ;
1691          T4_dVb = 2.0 / T2 * T0_dVb - 2.0 / T2 / T2 * T0 * T2_dVb ;
1692 
1693          Fn_SL_CP( T9 , T4 , 0 , DEPQFN_dlt, 2 , T0 )
1694          T9_dVg = T4_dVg * T0 ;
1695          T9_dVd = T4_dVd * T0 ;
1696          T9_dVb = T4_dVb * T0 ;
1697 
1698          T9 = T9 + small ;
1699          T3 = sqrt( T9 ) ;
1700          T3_dVg = 0.5 / T3 * T9_dVg  ;
1701          T3_dVd = 0.5 / T3 * T9_dVd  ;
1702          T3_dVb = 0.5 / T3 * T9_dVb  ;
1703 
1704          T10 = Vgp + here->HSM2_depvdsef1 + T2 * ( 1.0e0 - T3 ) ;
1705          T10 = T10 * here->HSM2_depvdsef2 ;
1706          T10_dVb = (Vgp_dVbs + T2_dVb * ( 1.0e0 - T3 ) - T2 * T3_dVb)
1707                    * here->HSM2_depvdsef2 ;
1708          T10_dVd = (Vgp_dVds + T2_dVd * ( 1.0e0 - T3 ) - T2 * T3_dVd)
1709                    * here->HSM2_depvdsef2 ;
1710          T10_dVg = (Vgp_dVgs + T2_dVg * ( 1.0e0 - T3 ) - T2 * T3_dVg)
1711                    * here->HSM2_depvdsef2 ;
1712 
1713          Fn_SL_CP( T10 , T10 , here->HSM2_depleak, 1.0, 4 , T0 )
1714          T10 = T10 + epsm10 ;
1715          T10_dVb *=T0 ;
1716          T10_dVd *=T0 ;
1717          T10_dVg *=T0 ;
1718 
1719          T1 = Vdsorg / T10 ;
1720          T2 = Fn_Pow( T1 , here->HSM2_ddlt - 1.0e0 ) ;
1721          T7 = T2 * T1 ;
1722          T0 = here->HSM2_ddlt * T2 / ( T10 * T10 ) ;
1723          T7_dVb = T0 * ( - Vdsorg * T10_dVb ) ;
1724          T7_dVd = T0 * ( T10 - Vdsorg * T10_dVd ) ;
1725          T7_dVg = T0 * ( - Vdsorg * T10_dVg ) ;
1726 
1727          T3 = 1.0 + T7 ;
1728          T4 = Fn_Pow( T3 , 1.0 / here->HSM2_ddlt - 1.0 ) ;
1729          T6 = T4 * T3 ;
1730          T0 = T4 / here->HSM2_ddlt ;
1731          T6_dVb = T0 * T7_dVb ;
1732          T6_dVd = T0 * T7_dVd ;
1733          T6_dVg = T0 * T7_dVg ;
1734 
1735          Vdseff0 = Vdsorg / T6 ;
1736          T0 = 1.0 / ( T6 * T6 ) ;
1737          Vdseff0_dVbs = - Vdsorg * T6_dVb * T0 ;
1738          Vdseff0_dVds = ( T6 - Vdsorg * T6_dVd ) * T0 ;
1739          Vdseff0_dVgs = - Vdsorg * T6_dVg * T0 ;
1740 
1741        } else {
1742 
1743          Vdseff0 = Vdsorg ;
1744          Vdseff0_dVgs = 0.0 ;
1745          Vdseff0_dVds = 1.0 ;
1746          Vdseff0_dVbs = 0.0 ;
1747 
1748        }
1749 
1750        T0 = here->HSM2_ninvd ;
1751 
1752        T4 = 1.0 + ( Pdsz ) * T0 ;
1753        T4_dVb = ( Pdsz_dVbs ) * T0 ;
1754        T4_dVd = ( Pdsz_dVds ) * T0 ;
1755        T4_dVg = ( Pdsz_dVgs ) * T0 ;
1756 
1757        Qiu = - Qn_res0 ;
1758        Qiu_dVgs = - Qn_res0_dVgs ;
1759        Qiu_dVds = - Qn_res0_dVds ;
1760        Qiu_dVbs = - Qn_res0_dVbs ;
1761 
1762        T5     = Qiu ;
1763        T5_dVb = Qiu_dVbs ;
1764        T5_dVd = Qiu_dVds ;
1765        T5_dVg = Qiu_dVgs ;
1766 
1767        T3     = T5 / T4 ;
1768        T3_dVb = ( - T4_dVb * T5 + T4 * T5_dVb ) / T4 / T4 ;
1769        T3_dVd = ( - T4_dVd * T5 + T4 * T5_dVd ) / T4 / T4 ;
1770        T3_dVg = ( - T4_dVg * T5 + T4 * T5_dVg ) / T4 / T4 ;
1771 
1772        Eeff = T3 ;
1773        Eeff_dVbs = T3_dVb ;
1774        Eeff_dVds = T3_dVd ;
1775        Eeff_dVgs = T3_dVg ;
1776 
1777        T5 = Fn_Pow( Eeff , model->HSM2_depmueph0 - 1.0e0 ) ;
1778        T8 = T5 * Eeff ;
1779 
1780        T9 = C_QE * C_m2cm_p2 ;
1781        Rns = Qiu / T9 ;
1782 
1783        T1 = 1.0e0 / ( here->HSM2_depmue0 + here->HSM2_depmue1 * Rns / 1.0e11 )
1784          + here->HSM2_depmphn0 * T8  ;
1785 
1786        Muun = 1.0e0 / T1 ;
1787 
1788        T1 = 1.0e0 / ( T1 * T1 ) ;
1789        T2 = here->HSM2_depmue0 + here->HSM2_depmue1 * Rns / 1.0e11 ;
1790        T2 = 1.0e0 / ( T2 * T2 ) ;
1791        T3 = here->HSM2_depmphn1 * T5 ;
1792        T5 = - 1.0e-11 * here->HSM2_depmue1 / C_QE * T2 / C_m2cm_p2 ;
1793        Muun_dVbs = - ( T5 * Qiu_dVbs
1794                     + Eeff_dVbs * T3 ) * T1 ;
1795        Muun_dVds = - ( T5 * Qiu_dVds
1796                     + Eeff_dVds * T3 ) * T1 ;
1797        Muun_dVgs = - ( T5 * Qiu_dVgs
1798                     + Eeff_dVgs * T3 ) * T1 ;
1799 
1800       /*  Change to MKS unit */
1801        Muun /=C_m2cm_p2 ;
1802        Muun_dVbs /=C_m2cm_p2 ;
1803        Muun_dVds /=C_m2cm_p2 ;
1804        Muun_dVgs /=C_m2cm_p2 ;
1805 
1806        Edri = Vdseff0 / Lch ;
1807        Edri_dVgs = Vdseff0_dVgs / Lch ;
1808        Edri_dVds = Vdseff0_dVds / Lch ;
1809        Edri_dVbs = Vdseff0_dVbs / Lch ;
1810 
1811        T1 = Muun * Edri / here->HSM2_depvmax ;
1812        T1_dVgs = (Muun_dVgs * Edri + Muun * Edri_dVgs) / here->HSM2_depvmax ;
1813        T1_dVds = (Muun_dVds * Edri + Muun * Edri_dVds) / here->HSM2_depvmax ;
1814        T1_dVbs = (Muun_dVbs * Edri + Muun * Edri_dVbs) / here->HSM2_depvmax ;
1815 
1816        T1 = T1 + small ;
1817        T2 = Fn_Pow(T1,model->HSM2_depbb) ;
1818        T2_dVgs = model->HSM2_depbb * T1_dVgs / T1 * T2 ;
1819        T2_dVds = model->HSM2_depbb * T1_dVds / T1 * T2 ;
1820        T2_dVbs = model->HSM2_depbb * T1_dVbs / T1 * T2 ;
1821 
1822        T3 = 1.0 + T2 ;
1823        T4 = Fn_Pow(T3,1.0 / model->HSM2_depbb) ;
1824        T4_dVgs = 1.0 / model->HSM2_depbb * T2_dVgs / T3 * T4 ;
1825        T4_dVds = 1.0 / model->HSM2_depbb * T2_dVds / T3 * T4 ;
1826        T4_dVbs = 1.0 / model->HSM2_depbb * T2_dVbs / T3 * T4 ;
1827 
1828        Mu_res = Muun / T4 ;
1829        Mu_res_dVgs = Muun_dVgs / T4 - Muun / T4 / T4 * T4_dVgs ;
1830        Mu_res_dVds = Muun_dVds / T4 - Muun / T4 / T4 * T4_dVds ;
1831        Mu_res_dVbs = Muun_dVbs / T4 - Muun / T4 / T4 * T4_dVbs ;
1832 
1833        Ids_res = here->HSM2_weff_nf * ( - Qn_res0) * Mu_res * Edri ;
1834        Ids_res_dVgs = here->HSM2_weff_nf * ( - Qn_res0_dVgs * Mu_res * Edri
1835                - Qn_res0 * Mu_res_dVgs * Edri - Qn_res0 * Mu_res * Edri_dVgs) ;
1836        Ids_res_dVds = here->HSM2_weff_nf * ( - Qn_res0_dVds * Mu_res * Edri
1837                - Qn_res0 * Mu_res_dVds * Edri - Qn_res0 * Mu_res * Edri_dVds) ;
1838        Ids_res_dVbs = here->HSM2_weff_nf * ( - Qn_res0_dVbs * Mu_res * Edri
1839                - Qn_res0 * Mu_res_dVbs * Edri - Qn_res0 * Mu_res * Edri_dVbs) ;
1840 
1841 
1842      //-----------------------------------------------------------*
1843      //*  back region universal mobility.  (CGS unit)
1844      //*-----------------//
1845 
1846        T0 = here->HSM2_ninvd ;
1847 
1848        T4 = 1.0 + ( Pdsz ) * T0 ;
1849        T4_dVb = ( Pdsz_dVbs ) * T0 ;
1850        T4_dVd = ( Pdsz_dVds ) * T0 ;
1851        T4_dVg = ( Pdsz_dVgs ) * T0 ;
1852 
1853        Qiu = - Qn_bac0 ;
1854        Qiu_dVgs = - Qn_bac0_dVgs ;
1855        Qiu_dVds = - Qn_bac0_dVds ;
1856        Qiu_dVbs = - Qn_bac0_dVbs ;
1857 
1858        T5     = Qiu ;
1859        T5_dVb = Qiu_dVbs ;
1860        T5_dVd = Qiu_dVds ;
1861        T5_dVg = Qiu_dVgs ;
1862 
1863        T3     = T5 / T4 ;
1864        T3_dVb = ( - T4_dVb * T5 + T4 * T5_dVb ) / T4 / T4 ;
1865        T3_dVd = ( - T4_dVd * T5 + T4 * T5_dVd ) / T4 / T4 ;
1866        T3_dVg = ( - T4_dVg * T5 + T4 * T5_dVg ) / T4 / T4 ;
1867 
1868        Eeff = T3 ;
1869        Eeff_dVbs = T3_dVb ;
1870        Eeff_dVds = T3_dVd ;
1871        Eeff_dVgs = T3_dVg ;
1872 
1873        T5 = Fn_Pow( Eeff , model->HSM2_depmueph0 - 1.0e0 ) ;
1874        T8 = T5 * Eeff ;
1875 
1876        T9 = C_QE * C_m2cm_p2 ;
1877        Rns = Qiu / T9 ;
1878 
1879        T1 = 1.0e0 / ( here->HSM2_depmueback0 + here->HSM2_depmueback1 * Rns / 1.0e11 )
1880          + here->HSM2_depmphn0 * T8  ;
1881 
1882        Muun = 1.0e0 / T1 ;
1883 
1884        T1 = 1.0e0 / ( T1 * T1 ) ;
1885        T2 = here->HSM2_depmueback0 + here->HSM2_depmueback1 * Rns / 1.0e11 ;
1886        T2 = 1.0e0 / ( T2 * T2 ) ;
1887        T3 = here->HSM2_depmphn1 * T5 ;
1888        T5 = - 1.0e-11 * here->HSM2_depmueback1 / C_QE * T2 / C_m2cm_p2 ;
1889        Muun_dVbs = - ( T5 * Qiu_dVbs
1890                     + Eeff_dVbs * T3 ) * T1 ;
1891        Muun_dVds = - ( T5 * Qiu_dVds
1892                     + Eeff_dVds * T3 ) * T1 ;
1893        Muun_dVgs = - ( T5 * Qiu_dVgs
1894                     + Eeff_dVgs * T3 ) * T1 ;
1895 
1896       /*  Change to MKS unit */
1897        Muun /=C_m2cm_p2 ;
1898        Muun_dVbs /=C_m2cm_p2 ;
1899        Muun_dVds /=C_m2cm_p2 ;
1900        Muun_dVgs /=C_m2cm_p2 ;
1901 
1902        Edri = Vdseff0 / Lch ;
1903        Edri_dVgs = Vdseff0_dVgs / Lch ;
1904        Edri_dVds = Vdseff0_dVds / Lch ;
1905        Edri_dVbs = Vdseff0_dVbs / Lch ;
1906 
1907        T1 = Muun * Edri / here->HSM2_depvmax ;
1908        T1_dVgs = (Muun_dVgs * Edri + Muun * Edri_dVgs) / here->HSM2_depvmax ;
1909        T1_dVds = (Muun_dVds * Edri + Muun * Edri_dVds) / here->HSM2_depvmax ;
1910        T1_dVbs = (Muun_dVbs * Edri + Muun * Edri_dVbs) / here->HSM2_depvmax ;
1911 
1912        T1 = T1 + small ;
1913        T2 = Fn_Pow(T1,model->HSM2_depbb) ;
1914        T2_dVgs = model->HSM2_depbb * T1_dVgs / T1 * T2 ;
1915        T2_dVds = model->HSM2_depbb * T1_dVds / T1 * T2 ;
1916        T2_dVbs = model->HSM2_depbb * T1_dVbs / T1 * T2 ;
1917 
1918        T3 = 1.0 + T2 ;
1919        T4 = Fn_Pow(T3,1.0 / model->HSM2_depbb) ;
1920        T4_dVgs = 1.0 / model->HSM2_depbb * T2_dVgs / T3 * T4 ;
1921        T4_dVds = 1.0 / model->HSM2_depbb * T2_dVds / T3 * T4 ;
1922        T4_dVbs = 1.0 / model->HSM2_depbb * T2_dVbs / T3 * T4 ;
1923 
1924 
1925        Mu_bac = Muun / T4 ;
1926        Mu_bac_dVgs = Muun_dVgs / T4 - Muun / T4 / T4 * T4_dVgs ;
1927        Mu_bac_dVds = Muun_dVds / T4 - Muun / T4 / T4 * T4_dVds ;
1928        Mu_bac_dVbs = Muun_dVbs / T4 - Muun / T4 / T4 * T4_dVbs ;
1929 
1930        Ids_bac = here->HSM2_weff_nf * ( - Qn_bac0) * Mu_bac * Edri ;
1931        Ids_bac_dVgs = here->HSM2_weff_nf * ( - Qn_bac0_dVgs * Mu_bac * Edri
1932                - Qn_bac0 * Mu_bac_dVgs * Edri - Qn_bac0 * Mu_bac * Edri_dVgs) ;
1933        Ids_bac_dVds = here->HSM2_weff_nf * ( - Qn_bac0_dVds * Mu_bac * Edri
1934                - Qn_bac0 * Mu_bac_dVds * Edri - Qn_bac0 * Mu_bac * Edri_dVds) ;
1935        Ids_bac_dVbs = here->HSM2_weff_nf * ( - Qn_bac0_dVbs * Mu_bac * Edri
1936                - Qn_bac0 * Mu_bac_dVbs * Edri - Qn_bac0 * Mu_bac * Edri_dVbs) ;
1937 
1938 
1939      /*-----------------------------------------------------------*
1940       * Ids: channel current.
1941       *-----------------*/
1942        betaWL = here->HSM2_weff_nf * beta_inv / Lch ;
1943        T1 = - betaWL / Lch ;
1944 
1945        Ids0 = betaWL * Idd * Mu + Ids_res + Ids_bac ;
1946 
1947        Ids0_dVgs = betaWL * ( Idd_dVgs * Mu + Idd * Mu_dVgs )
1948                    + Ids_res_dVgs + Ids_bac_dVgs ;
1949        Ids0_dVds = betaWL * ( Idd_dVds * Mu + Idd * Mu_dVds )
1950                    + Ids_res_dVds + Ids_bac_dVds ;
1951        Ids0_dVbs = betaWL * ( Idd_dVbs * Mu + Idd * Mu_dVbs )
1952                    + Ids_res_dVbs + Ids_bac_dVbs ;
1953 
1954 
1955        // Vdseff //
1956 
1957        Vds = Vdsorg;
1958 
1959       /*-----------------------------------------------------------*
1960        * Adding parasitic components to the channel current.
1961        *-----------------*/
1962        if( model->HSM2_ptl != 0 ){
1963          T1 =  0.5 * ( Vds - Pds ) ;
1964          Fn_SymAdd( T6 , T1 , 0.01 , T2 ) ;
1965          T2 = T2 * 0.5 ;
1966          T6_dVb = T2 * ( - Pds_dVbs ) ;
1967          T6_dVd = T2 * ( 1.0 - Pds_dVds ) ;
1968          T6_dVg = T2 * ( - Pds_dVgs ) ;
1969 
1970          T1     = 1.1 - ( phi_s0_DEP + T6 );
1971          T1_dVb = - ( phi_s0_DEP_dVbs + T6_dVb );
1972          T1_dVd = - ( phi_s0_DEP_dVds + T6_dVd );
1973          T1_dVg = - ( phi_s0_DEP_dVgs + T6_dVg );
1974 
1975          Fn_SZ( T2 , T1 , 0.05 , T0 ) ;
1976          T2 = T2 + small ;
1977          T2_dVb = T1_dVb * T0 ;
1978          T2_dVd = T1_dVd * T0 ;
1979          T2_dVg = T1_dVg * T0 ;
1980 
1981          T0 = beta * here->HSM2_ptl0 ;
1982          T3 = Cox * T0 ;
1983          T3_dVb = Cox_dVb * T0 ;
1984          T3_dVd = Cox_dVd * T0 ;
1985          T3_dVg = Cox_dVg * T0 ;
1986          T0 = pow( T2 , model->HSM2_ptp ) ;
1987          T9     = T3 * T0 ;
1988          T9_dVb = T3 * model->HSM2_ptp * T0 / T2 * T2_dVb + T3_dVb * T0 ;
1989          T9_dVd = T3 * model->HSM2_ptp * T0 / T2 * T2_dVd + T3_dVd * T0 ;
1990          T9_dVg = T3 * model->HSM2_ptp * T0 / T2 * T2_dVg + T3_dVg * T0 ;
1991 
1992 
1993          T4 = 1.0 + Vdsz * model->HSM2_pt2 ;
1994          T4_dVb = Vdsz_dVbs * model->HSM2_pt2 ;
1995          T4_dVd = Vdsz_dVds * model->HSM2_pt2 ;
1996          T4_dVg = 0.0 ;
1997 
1998          T0 = here->HSM2_pt40 ;
1999          T5 = phi_s0_DEP + T6 - Vbsz ;
2000          T5_dVb = phi_s0_DEP_dVbs + T6_dVb - Vbsz_dVbs ;
2001          T5_dVd = phi_s0_DEP_dVds + T6_dVd - Vbsz_dVds ;
2002          T5_dVg = phi_s0_DEP_dVgs + T6_dVg ;
2003          T4 = T4 + Vdsz * T0 * T5 ;
2004          T4_dVb = T4_dVb + Vdsz * T0 * T5_dVb + Vdsz_dVbs * T0 * T5 ;
2005          T4_dVd = T4_dVd + Vdsz * T0 * T5_dVd + Vdsz_dVds * T0 * T5 ;
2006          T4_dVg = T4_dVg + Vdsz * T0 * T5_dVg ;
2007          T6     = T9 * T4 ;
2008          T9_dVb = T9_dVb * T4 + T9 * T4_dVb ;
2009          T9_dVd = T9_dVd * T4 + T9 * T4_dVd ;
2010          T9_dVg = T9_dVg * T4 + T9 * T4_dVg ;
2011          T9     = T6 ;
2012 
2013        } else {
2014          T9 = 0.0 ;
2015          T9_dVb = 0.0 ;
2016          T9_dVd = 0.0 ;
2017          T9_dVg = 0.0 ;
2018        }
2019 
2020        if( model->HSM2_gdl != 0 ){
2021          T1 = beta * here->HSM2_gdl0 ;
2022          T2 = Cox * T1 ;
2023          T2_dVb = Cox_dVb * T1 ;
2024          T2_dVd = Cox_dVd * T1 ;
2025          T2_dVg = Cox_dVg * T1 ;
2026          T8     = T2 * Vdsz ;
2027          T8_dVb = T2_dVb * Vdsz + T2 * Vdsz_dVbs ;
2028          T8_dVd = T2_dVd * Vdsz + T2 * Vdsz_dVds ;
2029          T8_dVg = T2_dVg * Vdsz ;
2030        } else {
2031          T8 = 0.0 ;
2032          T8_dVb = 0.0 ;
2033          T8_dVd = 0.0 ;
2034          T8_dVg = 0.0 ;
2035        }
2036 
2037        if ( ( T9 + T8 ) > 0.0 ) {
2038          Idd1 = Pds * ( T9 + T8 ) ;
2039          Idd1_dVbs = Pds_dVbs * ( T9 + T8 ) + Pds * ( T9_dVb + T8_dVb ) ;
2040          Idd1_dVds = Pds_dVds * ( T9 + T8 ) + Pds * ( T9_dVd + T8_dVd ) ;
2041          Idd1_dVgs = Pds_dVgs * ( T9 + T8 ) + Pds * ( T9_dVg + T8_dVg ) ;
2042 
2043          Ids0 = Ids0 +betaWL * Idd1 * Mu ;
2044          T1 = betaWL * Idd1 ;
2045          T2 = Idd1 * Mu ;
2046          T3 = Mu * betaWL ;
2047          Ids0_dVbs +=T3 * Idd1_dVbs + T1 * Mu_dVbs + T2 * betaWL_dVbs ;
2048          Ids0_dVds +=T3 * Idd1_dVds + T1 * Mu_dVds + T2 * betaWL_dVds ;
2049          Ids0_dVgs +=T3 * Idd1_dVgs + T1 * Mu_dVgs + T2 * betaWL_dVgs ;
2050        }
2051 
2052        /* Ids calculation */
2053       if ( flg_rsrd == 2 ) {
2054          Rd  = here->HSM2_rd ;
2055          T0 = Rd * Ids0 ;
2056          T1 = Vds + small ;
2057          T2 = 1.0 / T1 ;
2058          T3 = 1.0 + T0 * T2 ;
2059          T3_dVb =   Rd * Ids0_dVbs             * T2 ;
2060          T3_dVd = ( Rd * Ids0_dVds * T1 - T0 ) * T2 * T2 ;
2061          T3_dVg =   Rd * Ids0_dVgs             * T2 ;
2062          T4 = 1.0 / T3 ;
2063          Ids = Ids0 * T4 ;
2064          T5 = T4 * T4 ;
2065          Ids_dVbs = ( Ids0_dVbs * T3 - Ids0 * T3_dVb ) * T5 ;
2066          Ids_dVds = ( Ids0_dVds * T3 - Ids0 * T3_dVd ) * T5 ;
2067          Ids_dVgs = ( Ids0_dVgs * T3 - Ids0 * T3_dVg ) * T5 ;
2068        } else {
2069          Ids = Ids0 ;
2070          Ids_dVbs = Ids0_dVbs ;
2071          Ids_dVds = Ids0_dVds ;
2072          Ids_dVgs = Ids0_dVgs ;
2073        }
2074 
2075       /* charge calculation */
2076       /*---------------------------------------------------*
2077        * Qbu : - Qb in unit area.
2078        *-----------------*/
2079        Qbu = - 0.5 * (Q_sub0 + Q_subL + Q_sub0_dep + Q_subL_dep ) ;
2080        Qbu_dVgs = - 0.5 * ( Q_sub0_dVgs + Q_subL_dVgs + Q_sub0_dep_dVgs + Q_subL_dep_dVgs ) ;
2081        Qbu_dVds = - 0.5 * ( Q_sub0_dVds + Q_subL_dVds + Q_sub0_dep_dVds + Q_subL_dep_dVds ) ;
2082        Qbu_dVbs = - 0.5 * ( Q_sub0_dVbs + Q_subL_dVbs + Q_sub0_dep_dVbs + Q_subL_dep_dVbs ) ;
2083 
2084        Qiu = - 0.5 * (Q_n0 + Q_nL + Q_s0_dep + Q_sL_dep + Q_b0_dep + Q_bL_dep) ;
2085        Qiu_dVgs = - 0.5 * ( Q_n0_dVgs + Q_nL_dVgs + Q_s0_dep_dVgs + Q_sL_dep_dVgs + Q_b0_dep_dVgs + Q_bL_dep_dVgs ) ;
2086        Qiu_dVds = - 0.5 * ( Q_n0_dVds + Q_nL_dVds + Q_s0_dep_dVds + Q_sL_dep_dVds + Q_b0_dep_dVds + Q_bL_dep_dVds ) ;
2087        Qiu_dVbs = - 0.5 * ( Q_n0_dVbs + Q_nL_dVbs + Q_s0_dep_dVbs + Q_sL_dep_dVbs + Q_b0_dep_dVbs + Q_bL_dep_dVbs ) ;
2088 
2089        Qdrat = 0.5;
2090        Qdrat_dVgs = 0.0 ;
2091        Qdrat_dVds = 0.0 ;
2092        Qdrat_dVbs = 0.0 ;
2093 
2094       /*-------------------------------------------------*
2095        * set flg_noqi , Qn0 , Ey for noise calc.
2096        *-----------------*/
2097        Qiu_noi = - 0.5 * (Q_n0 + Q_nL ) ;
2098        Qn0 = - Q_n0 ;
2099        Qn0_dVgs = - Q_n0_dVgs ;
2100        Qn0_dVds = - Q_n0_dVds ;
2101        Qn0_dVbs = - Q_n0_dVbs ;
2102 
2103        Ey = Ey_suf ;
2104        Ey_dVbs = Ey_suf_dVbs ;
2105        Ey_dVds = Ey_suf_dVds ;
2106        Ey_dVgs = Ey_suf_dVgs ;
2107 
2108        if( Qn0 < small ) { flg_noqi = 1; }
2109 
2110 } /* End of hsmhveval_dep */
2111 
2112 
2113