mkfact.f
189 lines
| 3.9 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | ||||
| subroutine mkfact(tm,h,bfm,thb,bki,nhts) | ||||
| c | ||||
| c create table of magnetic field components | ||||
| c | ||||
| real h(nhts),bfm(nhts),thb(nhts),bki(nhts) | ||||
| data pio180/0.017453292/,pi/3.14159265/ | ||||
| C rate of rotation produced by n=1 cm-3, b=1 gauss | ||||
| data rate/1.8978873e-6/ | ||||
| C Jicamarca position | ||||
| data gdlat/-11.951/,elon/-76.8667/,gdalt/0./,re/6280.0/ | ||||
| C 3.0, 4.5, 6.0 and 3.38N degree positions | ||||
| c data decd/-12.88/,ham/-4.617/ | ||||
| c data decd/-13.1/,ham/-5.0/ ! made up | ||||
| c data decd/-14.56/,ham/-5.33/ | ||||
| c data decd/-16.24/,ham/-6.03/ | ||||
| data decd/-9.52/,ham/-3.20/ | ||||
| c======= | ||||
| c | ||||
| c convert from geodetic to geocentric coordinates | ||||
| c | ||||
| call convrt(1,gdlat,gdalt,gclat,re) | ||||
| c | ||||
| c convert from degrees to radians | ||||
| c | ||||
| gclat=gclat*pio180 | ||||
| elong=elon*pio180 | ||||
| c | ||||
| c create components of vector from center of earth | ||||
| c to antenna center | ||||
| ca=cos(elong) | ||||
| cd=cos(gclat) | ||||
| sa=sin(elong) | ||||
| sd=sin(gclat) | ||||
| rx=re*ca*cd | ||||
| ry=re*sa*cd | ||||
| rz=re*sd | ||||
| c | ||||
| c create antenna normal vector (unit vector in direction of k) | ||||
| c | ||||
| dec=decd*pio180 | ||||
| ha=ham*pi/720.+elong | ||||
| ch=cos(ha) | ||||
| sh=sin(ha) | ||||
| cdec=cos(dec) | ||||
| sdec=sin(dec) | ||||
| ax=ch*cdec | ||||
| ay=sh*cdec | ||||
| az=sdec | ||||
| do i=1,nhts | ||||
| c | ||||
| c create vector from center of earth to scattering volume | ||||
| c | ||||
| x=rx+ax*h(i) | ||||
| y=ry+ay*h(i) | ||||
| z=rz+az*h(i) | ||||
| x2y2=x*x+y*y | ||||
| r=sqrt(x2y2+z*z) | ||||
| x2y2=sqrt(x2y2) | ||||
| c | ||||
| c convert to spherical coordinates | ||||
| c | ||||
| st=x2y2/r | ||||
| ct=z/r | ||||
| sp=y/x2y2 | ||||
| cp=x/x2y2 | ||||
| c call milmag(tm,r,st,ct,sp,cp,br,bt,bp,b) | ||||
| theta=atan2(st,ct) | ||||
| phi=atan2(sp,cp) | ||||
| call geobfield(tm,r,theta,phi,br,bt,bp,b) | ||||
| c | ||||
| c convert to cartesian coordinates | ||||
| c | ||||
| bx=br*st*cp+bt*ct*cp-bp*sp | ||||
| by=br*st*sp+bt*ct*sp+bp*cp | ||||
| bz=br*ct-bt*st | ||||
| bk=(ax*bx+ay*by+az*bz) | ||||
| thb(i)=acos(abs(bk/b))*57.29577951 | ||||
| bfm(i)=b | ||||
| bki(i)=1.0/bk | ||||
| end do | ||||
| return | ||||
| end | ||||
| subroutine gausfit(drift,hist,vz,m) | ||||
| c | ||||
| c fit procedure for estimating velocities | ||||
| c | ||||
| real drift(m),hist(m),vz | ||||
| real x(1000),y(1000) | ||||
| parameter(n=6) | ||||
| real a(n),tol,wa(m*n+5*n+m),wa2(m*n+5*n+m),fvec(m) | ||||
| integer info,iwa(n),lwa | ||||
| external gauss | ||||
| common/gf/x,y | ||||
| data vlast/0.0/ | ||||
| lwa=m*n+5*n+m | ||||
| tol=1.0e-6 | ||||
| info=0 | ||||
| amax=0.0 | ||||
| do i=1,m | ||||
| x(i)=drift(i) | ||||
| y(i)=hist(i) | ||||
| amax=max(amax,hist(i)) | ||||
| end do | ||||
| a(1)=amax | ||||
| a(2)=-20.0 | ||||
| a(3)=20.0 | ||||
| a(4)=amax/4.0 | ||||
| a(5)=20.0 | ||||
| a(6)=20.0 | ||||
| c write(*,*) m,n,a,tol,info | ||||
| call lmdif1(gauss,m,n,a,fvec,tol,info,iwa,wa,lwa) | ||||
| call fdjac2(gauss,m,n,a,fvec,wa,m,iflag,0.0e0,wa2) | ||||
| if(abs(a(1)).gt.abs(a(4))) then | ||||
| vz=a(2) | ||||
| i=2 | ||||
| else | ||||
| vz=a(5) | ||||
| i=5 | ||||
| end if | ||||
| err=0.0 | ||||
| do j=1,m | ||||
| err=err+(wa(j+(i-1)*m))**2 | ||||
| end do | ||||
| if(err.gt.0.0) err=sqrt(err**-1) | ||||
| c | ||||
| c try fitting single distribution | ||||
| c | ||||
| a(1)=amax | ||||
| a(2)=-20.0 | ||||
| a(3)=20.0 | ||||
| call lmdif1(gauss,m,n/2,a,fvec,tol,info,iwa,wa,lwa) | ||||
| call fdjac2(gauss,m,n/2,a,fvec,wa,m,iflag,0.0e0,wa2) | ||||
| vz1=a(2) | ||||
| i=2 | ||||
| err1=0.0 | ||||
| do j=1,m | ||||
| err1=err1+(wa(j+(i-1)*m))**2 | ||||
| end do | ||||
| if(err1.gt.0.0) err1=sqrt(err**-1) | ||||
| if(err1.lt.err.or.abs(vz-vz1).lt.10.0) vz=vz1 | ||||
| vz=vz1 | ||||
| vlast=vz | ||||
| return | ||||
| end | ||||
| subroutine gauss(m,n,a,fvec,iflag) | ||||
| c | ||||
| integer m,n,iflac | ||||
| real fvec(m),a(n),b | ||||
| real x(1000),y(1000) | ||||
| common/gf/x,y | ||||
| c | ||||
| do i=1,m | ||||
| sig=sqrt(max(1.0,y(i))) | ||||
| b=abs(a(1))*exp(-(a(2)-x(i))**2/(2.0*a(3)**2)) | ||||
| if(n.gt.3) then | ||||
| c=abs(a(4))*exp(-(a(5)-x(i))**2/(2.0*a(6)**2)) | ||||
| else | ||||
| c=0.0 | ||||
| end if | ||||
| fvec(i)=(b+c-y(i))**2/sig**2 | ||||
| end do | ||||
| return | ||||
| end | ||||
