Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_BNS_peos_lecc

00001 subroutine update_parameter_BNS_peos_lecc_spin(convf)
00002   use phys_constant, only  : long
00003   use grid_parameter, only : nrg, ntg, npg, nrf_deform, nrf, ntgxy, npgxzm
00004   use def_metric, only : psi, alph, bvxd, bvyd, bvzd, tfkij, tfkijkij, alps
00005   use def_matter_parameter
00006   use def_velocity_rot
00007   use def_velocity_potential
00008   use def_matter
00009   use coordinate_grav_r, only : rg
00010   use trigonometry_grav_theta, only : sinthg, costhg
00011   use trigonometry_grav_phi, only   : sinphig, cosphig
00012   use def_vector_phi, only : vec_phig
00013   use interface_minv
00014   use interface_flgrad_4th_gridpoint
00015   implicit none
00016   real(long), intent(in) :: convf
00017   real(long) ::    sgg(1:4), gg(1:4,1:4)
00018   real(long) ::    ain, amid, aout, pin, pmid, pout
00019 ! ain is log(\alpha^((1/R_i)^2)) (in). R_in is in fact the previous R.
00020   real(long) :: a2ppin,  a2p4pin,  pm4in,  ppin,  p4pin,  bxin,  byin,  bzin
00021   real(long) :: a2ppmid, a2p4pmid, pm4mid, ppmid, p4pmid, bxmid, bymid, bzmid  
00022   real(long) :: a2ppout, a2p4pout, pm4out, ppout, p4pout, bxout, byout, bzout 
00023   real(long) :: dxvepin,  dyvepin,  dzvepin,  wxin,  wyin,  wzin
00024   real(long) :: dxvepmid, dyvepmid, dzvepmid, wxmid, wymid, wzmid
00025   real(long) :: dxvepout, dyvepout, dzvepout, wxout, wyout, wzout
00026   real(long) :: ovin, ovmid, ovout, lamin, lammid, lamout, lam2in, lam2mid, lam2out
00027   real(long) :: vphidvepin, vphidvepmid, vphidvepout
00028   real(long) :: dvep2in,  wdvepin,  w2in,  term1in,  term2in,  term3in,  term4in,  term5in,  term6in
00029   real(long) :: dvep2mid, wdvepmid, w2mid, term1mid, term2mid, term3mid, term4mid, term5mid, term6mid
00030   real(long) :: dvep2out, wdvepout, w2out, term1out, term2out, term3out, term4out, term5out, term6out
00031   real(long) :: vphiyin, vphiymid, vphiyout
00032   real(long) :: loghc
00033   real(long) :: radiold, berold, omeold
00034   real(long) :: ddradi, ddber, ddome
00035   real(long) :: facfac, numero, error, hc, pre, rho, ene
00036   integer    :: nin
00037   integer    :: irin, itin, ipin, irmid, itmid, ipmid, irout, itout, ipout
00038 !##############################################
00039   real(long) ::    ut, hh, emdtest
00040 !##############################################
00041   nin = nrf_deform
00042 !
00043   sgg = 0.0d0
00044   gg  = 0.0d0
00045 !
00046   irin = nrf;   itin = ntgxy;   ipin = 0;
00047   irmid= 0  ;   itmid= 0    ;   ipmid= 0;
00048   irout= nrf;   itout= ntgxy;   ipout= npgxzm;
00049 
00050   call peos_q2hprho(emdc, hc, pre, rho, ene)
00051   loghc = dlog(hc)
00052 
00053   vphiyin  = vec_phig(irin,  itin,  ipin,  2)
00054   vphiymid = vec_phig(irmid, itmid, ipmid, 2)
00055   vphiyout = vec_phig(irout, itout, ipout, 2)
00056 
00057   omeold = ome
00058   berold = ber
00059   radiold = radi
00060 
00061 !  write(6,'(1p,6e12.4)') ome, ber, radi, vphiyin, vphiymid, vphiyout
00062 
00063   numero = 1
00064 ! ################################################################
00065   ain     = dlog(alph(irin,itin,ipin)**(1.d0/radi**2))
00066   pin     = dlog( psi(irin,itin,ipin)**(1.d0/radi**2))
00067   a2ppin  = ( alph(irin,itin,ipin)**2 * psi(irin,itin,ipin)**(confpow)               )**(1.d0/radi**2)
00068   a2p4pin = ( alph(irin,itin,ipin)**2 * psi(irin,itin,ipin)**(4.0d0 + 2.0d0*confpow) )**(1.d0/radi**2)
00069   pm4in   = ( psi(irin,itin,ipin)**(-4)                                              )**(1.d0/radi**2)
00070   ppin    = ( psi(irin,itin,ipin)**(confpow)                                         )**(1.d0/radi**2)
00071   p4pin   = ( psi(irin,itin,ipin)**(4.0d0 + 2.0d0*confpow)                           )**(1.d0/radi**2)
00072   bxin    = bvxd(irin,itin,ipin)
00073   byin    = bvyd(irin,itin,ipin)
00074   bzin    = bvzd(irin,itin,ipin)
00075 !
00076 !  write(6,*)  ain, pin, confpow, a2ppin
00077 
00078   amid     = dlog(alph(irmid,itmid,ipmid)**(1.d0/radi**2))
00079   pmid     = dlog( psi(irmid,itmid,ipmid)**(1.d0/radi**2))
00080   a2ppmid  = ( alph(irmid,itmid,ipmid)**2 * psi(irmid,itmid,ipmid)**(confpow)               )**(1.d0/radi**2)
00081   a2p4pmid = ( alph(irmid,itmid,ipmid)**2 * psi(irmid,itmid,ipmid)**(4.0d0 + 2.0d0*confpow) )**(1.d0/radi**2)
00082   pm4mid   = ( psi(irmid,itmid,ipmid)**(-4)                                                 )**(1.d0/radi**2)
00083   ppmid    = ( psi(irmid,itmid,ipmid)**(confpow)                                            )**(1.d0/radi**2)
00084   p4pmid   = ( psi(irmid,itmid,ipmid)**(4.0d0 + 2.0d0*confpow)                              )**(1.d0/radi**2)
00085   bxmid    = bvxd(irmid,itmid,ipmid)
00086   bymid    = bvyd(irmid,itmid,ipmid)
00087   bzmid    = bvzd(irmid,itmid,ipmid)
00088 !
00089   aout     = dlog(alph(irout,itout,ipout)**(1.d0/radi**2))
00090   pout     = dlog( psi(irout,itout,ipout)**(1.d0/radi**2))
00091   a2ppout  = ( alph(irout,itout,ipout)**2 * psi(irout,itout,ipout)**(confpow)               )**(1.d0/radi**2)
00092   a2p4pout = ( alph(irout,itout,ipout)**2 * psi(irout,itout,ipout)**(4.0d0 + 2.0d0*confpow) )**(1.d0/radi**2)
00093   pm4out   = ( psi(irout,itout,ipout)**(-4)                                                 )**(1.d0/radi**2)
00094   ppout    = ( psi(irout,itout,ipout)**(confpow)                                            )**(1.d0/radi**2)
00095   p4pout   = ( psi(irout,itout,ipout)**(4.0d0 + 2.0d0*confpow)                              )**(1.d0/radi**2)
00096   bxout    = bvxd(irout,itout,ipout)
00097   byout    = bvyd(irout,itout,ipout)
00098   bzout    = bvzd(irout,itout,ipout)
00099 
00100 ! write (6,*) "test values"
00101 ! write (6,'(1p,3e12.4)') ome, ber, radi
00102 ! write (6,'(1p,3e12.4)') psi(irin,itin,ipin) ,psi(irmid,itmid,ipmid) ,psi(irout,itout,ipout)
00103 ! write (6,'(1p,3e12.4)') alph(irin,itin,ipin),alph(irmid,itmid,ipmid),alph(irout,itout,ipout)
00104 ! write (6,'(1p,3e12.4)') error
00105 ! write (6,'(1p,3e12.4)') pain, pamid, paout
00106 !  write(6,'(a14,f10.7,a19,f10.7,a12,f10.7,a7,f10.7,a6)') &
00107 !  &   "ff=-Log[ber]+(",ain,")*R*R+0.5*Log[1.0-(",pain,")^(2*R*R)*((",bin,")+ome*(",vphiyin,"))^2];"
00108 !  write(6,'(a3,f10.7,a11,f10.7,a19,f10.7,a12,f10.7,a7,f10.7,a6)') "gg=",loghc,    &
00109 !  &   "-Log[ber]+(",amid,")*R*R+0.5*Log[1.0-(",pamid,")^(2*R*R)*((",bmid,")+ome*(",vphiymid,"))^2];"
00110 !  write(6,'(a14,f10.7,a19,f10.7,a12,f10.7,a7,f10.7,a6)') &
00111 !  &   "hh=-Log[ber]+(",aout,")*R*R+0.5*Log[1.0-(",paout,")^(2*R*R)*((",bout,")+ome*(",vphiyout,"))^2];"
00112 !  write(6,'(a7,f10.7,a7,f10.7,a5,f10.7,a2)') ",{{ome,",ome,"},{ber,",ber,"},{R,",radi,"}}"
00113 ! ################################################################
00114   do
00115     dxvepin = vepxf(irin,itin,ipin)
00116     dyvepin = vepyf(irin,itin,ipin)
00117     dzvepin = vepzf(irin,itin,ipin)
00118    
00119     wxin    = wxspf(irin,itin,ipin)
00120     wyin    = wyspf(irin,itin,ipin)
00121     wzin    = wzspf(irin,itin,ipin)
00122 
00123     dxvepmid = vepxf(irmid,itmid,ipmid)
00124     dyvepmid = vepyf(irmid,itmid,ipmid)
00125     dzvepmid = vepzf(irmid,itmid,ipmid)
00126 
00127     wxmid    = wxspf(irmid,itmid,ipmid)
00128     wymid    = wyspf(irmid,itmid,ipmid)
00129     wzmid    = wzspf(irmid,itmid,ipmid)
00130 
00131     dxvepout = vepxf(irout,itout,ipout)
00132     dyvepout = vepyf(irout,itout,ipout)
00133     dzvepout = vepzf(irout,itout,ipout)
00134 
00135     wxout    = wxspf(irout,itout,ipout)
00136     wyout    = wyspf(irout,itout,ipout)
00137     wzout    = wzspf(irout,itout,ipout)
00138 
00139     ovin  = byin  + ome*vphiyin
00140     ovmid = bymid + ome*vphiymid
00141     ovout = byout + ome*vphiyout
00142 
00143 !    lamin  = ber + (bxin +0)*dxvepin  +  ovin*dyvepin  + (bzin +0)*dzvepin
00144 !    lammid = ber + (bxmid+0)*dxvepmid + ovmid*dyvepmid + (bzmid+0)*dzvepmid
00145 !    lamout = ber + (bxout+0)*dxvepout + ovout*dyvepout + (bzout+0)*dzvepout
00146 
00147     lamin  = ber + ovin*dyvepin 
00148     lammid = ber + ovmid*dyvepmid 
00149     lamout = ber + ovout*dyvepout
00150 
00151     lam2in  = lamin*lamin
00152     lam2mid = lammid*lammid
00153     lam2out = lamout*lamout
00154  
00155     vphidvepin  =  vphiyin*dyvepin
00156     vphidvepmid = vphiymid*dyvepmid
00157     vphidvepout = vphiyout*dyvepout
00158 
00159 !    write(6,*)  vphiyin, dxvepin, dyvepin, dzvepin
00160 
00161     dvep2in = (dxvepin**2 + dyvepin**2 + dzvepin**2)*pm4in**(radi**2)
00162 
00163     wdvepin = (wxin*dxvepin + wyin*dyvepin + wzin*dzvepin)*ppin**(radi**2)
00164 
00165     w2in    = (wxin*wxin + wyin*wyin + wzin*wzin)*p4pin**(radi**2)
00166 
00167     term1in = (wxin*dxvepin + wyin*dyvepin + wzin*dzvepin) * a2ppin**(radi**2) + &
00168             & (wxin*wxin + wyin*wyin + wzin*wzin) * a2p4pin**(radi**2)
00169 
00170     term2in = dvep2in + 2.0d0*wdvepin + w2in
00171 
00172     term3in = 0.25d0 + term1in/lam2in
00173 
00174     term4in = 1.0d0 + term2in/1.0d0
00175 
00176     dvep2mid = (dxvepmid**2 + dyvepmid**2 + dzvepmid**2)*pm4mid**(radi**2)
00177 
00178     wdvepmid = (wxmid*dxvepmid + wymid*dyvepmid + wzmid*dzvepmid)*ppmid**(radi**2)
00179 
00180     w2mid    = (wxmid*wxmid + wymid*wymid + wzmid*wzmid)*p4pmid**(radi**2)
00181 
00182     term1mid = (wxmid*dxvepmid + wymid*dyvepmid + wzmid*dzvepmid) * a2ppmid**(radi**2) + &
00183             & (wxmid*wxmid + wymid*wymid + wzmid*wzmid) * a2p4pmid**(radi**2)
00184 
00185     term2mid = dvep2mid + 2.0d0*wdvepmid + w2mid
00186 
00187     term3mid = 0.25d0 + term1mid/lam2mid
00188 
00189     term4mid = 1.0d0 + term2mid/hc/hc
00190 
00191     dvep2out = (dxvepout**2 + dyvepout**2 + dzvepout**2)*pm4out**(radi**2)
00192 
00193     wdvepout = (wxout*dxvepout + wyout*dyvepout + wzout*dzvepout)*ppout**(radi**2)
00194 
00195     w2out    = (wxout*wxout + wyout*wyout + wzout*wzout)*p4pout**(radi**2)
00196 
00197     term1out = (wxout*dxvepout + wyout*dyvepout + wzout*dzvepout) * a2ppout**(radi**2) + &
00198             & (wxout*wxout + wyout*wyout + wzout*wzout) * a2p4pout**(radi**2)
00199 
00200     term2out = dvep2out + 2.0d0*wdvepout + w2out
00201 
00202     term3out = 0.25d0 + term1out/lam2out
00203 
00204     term4out = 1.0d0 + term2out/1.0d0
00205 !
00206     if(lamin<0 .or. lamout<0 .or. lammid<0) then
00207       write(6,*)  "Quantities lam inside log negative. Exiting..."
00208       stop
00209     end if
00210     if(term3in<0 .or. term3out<0 .or. term3mid<0) then
00211       write(6,*)  "Quantities inside sqrt negative. Exiting..."
00212       stop
00213     end if
00214     if(term4in<0 .or. term4out<0 .or. term4mid<0) then
00215       write(6,*)  "Quantities inside log negative. Exiting..."
00216       stop
00217     else
00218       sgg(1) = -( -dlog(lamin)  + radi**2*ain  - dlog(0.5d0+sqrt(term3in))  + &
00219              &     0.5d0*dlog(term4in))
00220 
00221       sgg(2) = -( -dlog(lamout) + radi**2*aout - dlog(0.5d0+sqrt(term3out)) + &
00222              &     0.5d0*dlog(term4out))
00223 
00224       sgg(3) = -( -dlog(lammid) + radi**2*amid - dlog(0.5d0+sqrt(term3mid)) + &
00225                    0.5d0*dlog(term4mid) + loghc)
00226     end if
00227 !
00228     term5in  =  a2ppin**(radi**2)*2.0d0*radi*dlog(a2ppin) *         &
00229              &  (wxin*dxvepin + wyin*dyvepin + wzin*dzvepin) +      &
00230              &  a2p4pin**(radi**2)*2.0d0*radi*dlog(a2p4pin)*        &
00231              &  (wxin*wxin + wyin*wyin + wzin*wzin)
00232 
00233     term6in  = dvep2in*2.0d0*radi*dlog(pm4in) + 2.0d0*wdvepin*2.0d0*radi*dlog(ppin) + &
00234              & w2in*2.0d0*radi*dlog(p4pin)
00235 
00236     term5mid =  a2ppmid**(radi**2)*2.0d0*radi*dlog(a2ppmid) *             &
00237              &  (wxmid*dxvepmid + wymid*dyvepmid + wzmid*dzvepmid) +      &
00238              &  a2p4pmid**(radi**2)*2.0d0*radi*dlog(a2p4pmid)*             &
00239              &  (wxmid*wxmid + wymid*wymid + wzmid*wzmid)
00240 
00241     term6mid = dvep2mid*2.0d0*radi*dlog(pm4mid) + 2.0d0*wdvepmid*2.0d0*radi*dlog(ppmid) + &
00242              & w2mid*2.0d0*radi*dlog(p4pmid)
00243 
00244     term5out =  a2ppout**(radi**2)*2.0d0*radi*dlog(a2ppout) *               &
00245              &  (wxout*dxvepout + wyout*dyvepout + wzout*dzvepout) +         &
00246              &  a2p4pout**(radi**2)*2.0d0*radi*dlog(a2p4pout)*               &
00247                 (wxout*wxout + wyout*wyout + wzout*wzout)
00248 
00249     term6out = dvep2out*2.0d0*radi*dlog(pm4out) + 2.0d0*wdvepout*2.0d0*radi*dlog(ppout) + &
00250              & w2out*2.0d0*radi*dlog(p4pout)
00251 
00252 !  ain is ln(\alpha)
00253     gg(1,1) = vphidvepin*(-1.0d0 +      &
00254             & term1in/(0.5d0+sqrt(term3in))/sqrt(term3in)/lam2in )/lamin
00255 
00256     gg(1,2) = (-1.0d0 + term1in/(0.5d0+sqrt(term3in))/sqrt(term3in)/lam2in)/lamin
00257 
00258     gg(1,3) = 2.0d0*radi*ain - term5in/(0.5d0+sqrt(term3in))/sqrt(term3in)/(2.0d0*lam2in) + &
00259             & term6in/term4in/(2.0d0*1.0d0)
00260 
00261     gg(2,1) = vphidvepout*(-1.0d0 +     &
00262             & term1out/(0.5d0+sqrt(term3out))/sqrt(term3out)/lam2out )/lamout
00263 
00264     gg(2,2) = (-1.0d0 + term1out/(0.5d0+sqrt(term3out))/sqrt(term3out)/lam2out)/lamout
00265 
00266     gg(2,3) = 2.0d0*radi*aout - term5out/(0.5d0+sqrt(term3out))/sqrt(term3out)/(2.0d0*lam2out) + &
00267             & term6out/term4out/(2.0d0*1.0d0)
00268 
00269     gg(3,1) = vphidvepmid*(-1.0d0 +    &
00270             & term1mid/(0.5d0+sqrt(term3mid))/sqrt(term3mid)/lam2mid )/lammid
00271 
00272     gg(3,2) = (-1.0d0 + term1mid/(0.5d0+sqrt(term3mid))/sqrt(term3mid)/lam2mid)/lammid
00273 
00274     gg(3,3) = 2.0d0*radi*amid - term5mid/(0.5d0+sqrt(term3mid))/sqrt(term3mid)/(2.0d0*lam2mid) + &
00275             & term6mid/term4mid/(2.0d0*hc*hc)
00276 !
00277 !    test
00278 !if (numero.le.20) then
00279 ! write (6,'(1p,8e12.4)') sgg(1), sgg(2), sgg(3), pamid, radi, ov2mid, termmid
00280 ! write(6,*) "--------------------------------------------------------"
00281 ! write (6,'(1p,4f20.12)') gg(1,1), gg(1,2), gg(1,3), sgg(1)
00282 ! write (6,'(1p,4f20.12)') gg(2,1), gg(2,2), gg(2,3), sgg(2)
00283 ! write (6,'(1p,4f20.12)') gg(3,1), gg(3,2), gg(3,3), sgg(3)
00284 ! write(6,*) "--------------------------------------------------------"
00285 !end if
00286 !stop
00287 !!pause
00288 !    test
00289     call minv(gg,sgg,3,4)
00290 
00291 !    write (6,'(1p,4e12.4)') numero, sgg(1), sgg(2), sgg(3)
00292 !    write(6,*) "Old:", numero, ome, ber, radi
00293 
00294     ddome = sgg(1)
00295     ddber = sgg(2)
00296     ddradi = sgg(3)
00297     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00298     ome = ome + ddome*facfac
00299     ber = ber + ddber*facfac
00300     radi = radi + ddradi*facfac
00301 
00302 !    write(6,*) "New:", ome, ber, radi
00303     error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00304 !    numero = numero + 1
00305 !    if (numero<1000) then
00306 !      write(6,*)' numero = ', numero, '   error =',error
00307 !    end if
00308 !      if (iter<10.and.numero>10) exit
00309     if (numero>1010) exit
00310 !    if (dabs(error)<1.0d-08) exit
00311     if (dabs(error)<1.0d-12) then
00312 !      write(6,*) "New:", numero, ome, ber, radi
00313       exit
00314     end if
00315     numero = numero + 1
00316   end do
00317 !    test
00318 ! write (6,*) "test values"
00319 ! write (6,*) sgg(1), sgg(2), sgg(3)
00320 ! write (6,*) ome
00321 !    test
00322 !      write(6,*)' numero = ', numero, '  error = ',error
00323 !
00324   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00325   if (radi < 0.d0) radi = - radi
00326   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00327   if (ome < 0.d0) ome = - ome
00328   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00329   if (ber < 0.d0) ber = - ber
00330   ome = convf*ome + (1.d0-convf)*omeold
00331   ber = convf*ber + (1.d0-convf)*berold
00332   radi = convf*radi + (1.d0-convf)*radiold
00333 !
00334 !
00335 ! --  Improving alpha and psi.
00336 !
00337   alph = alph**((radi/radiold)**2)
00338   psi = psi**((radi/radiold)**2)
00339 
00340   alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00341 !##############################################
00342 !
00343 !  stop
00344 end subroutine update_parameter_BNS_peos_lecc_spin
00345
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_BNS_peos_lecc_spin.f90
