/[PAMELA software]/DarthVader/OrbitalInfo/src/igrf_sub.for
ViewVC logotype

Diff of /DarthVader/OrbitalInfo/src/igrf_sub.for

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.3 by mocchiut, Tue Aug 4 14:01:27 2009 UTC revision 1.6 by mocchiut, Tue May 15 14:31:14 2012 UTC
# Line 2  Line 2 
2       &          xl,icode,dip,dec)       &          xl,icode,dip,dec)
3  c----------------------------------------------------------------  c----------------------------------------------------------------
4  c   INPUT:  c   INPUT:
5  c       xlat    geodatic latitude in degrees  c     xlat     geodatic latitude in degrees
6  c       xlong   geodatic longitude in degrees  c     xlong    geodatic longitude in degrees
7  c       year    decimal year (year+month/12.0-0.5 or year+day-of-year/365  c     year     decimal year (year+month/12.0-0.5 or year+day-of-year/365
8  c               or 366 if leap year)  c          or 366 if leap year)
9  c       height  height in km  c     height     height in km
10  c   OUTPUT:  c   OUTPUT:
11  c       xl      L value  c     xl     L value
12  c       icode   =1  L is correct; =2  L is not correct;  c     icode     =1  L is correct; =2  L is not correct;
13  c               =3  an approximation is used  c          =3  an approximation is used
14  c       dip     geomagnetic inclination in degrees  c     dip     geomagnetic inclination in degrees
15  c       dec     geomagnetic declination in degress  c     dec     geomagnetic declination in degress
16  c----------------------------------------------------------------  c----------------------------------------------------------------
17    
18        REAL              LATI,LONGI        REAL              LATI,LONGI
# Line 20  c--------------------------------------- Line 20  c---------------------------------------
20        SAVE /GENER/        SAVE /GENER/
21  C  C
22        CALL INITIZE        CALL INITIZE
23          ibbb=0        ibbb=0
24        ALOG2=ALOG(2.)        ALOG2=ALOG(2.)
25        ISTART=1        ISTART=1
26          lati=xlat          lati=xlat
# Line 32  c Line 32  c
32          CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)          CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)
33          CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)          CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)
34          DIP=ASIN(BDOWN/BABS)/UMR          DIP=ASIN(BDOWN/BABS)/UMR
35          DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR        DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR
36        RETURN        RETURN
37        END        END
38  c  c
# Line 191  C                          APPROXIMATION Line 191  C                          APPROXIMATION
191  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS
192  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
193        DIMENSION         V(3),U(3,3),P(8,100),SP(3)        DIMENSION         V(3),U(3,3),P(8,100),SP(3)
194        COMMON            X(3),H(144)        COMMON            X(3),H(196)
195        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
196        SAVE /FIDB0/        SAVE /FIDB0/
197        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 342  C-- The minimal allowable value of FI wa Line 342  C-- The minimal allowable value of FI wa
342  C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.  C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.
343  C-- D. Bilitza, Nov 87.  C-- D. Bilitza, Nov 87.
344  C  C
345  11      FI=0.5*ABS(FI)/SQRT(B0)+1E-12                        11         FI=0.5*ABS(FI)/SQRT(B0)+1E-12                      
346  C  C
347  C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.    C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
348  C  C
349  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
350  C  C
351        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
352          arg1=alog(FI)        arg1=alog(FI)
353          arg2=alog(DIMOB0)        arg2=alog(DIMOB0)
354  c       arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
355  c       if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
356        XX=3*arg1-arg2        XX=3*arg1-arg2
357        IF(XX.GT.23.0) GOTO 776          IF(XX.GT.23.0) GOTO 776  
358        IF(XX.GT.11.7) GOTO 775          IF(XX.GT.11.7) GOTO 775  
# Line 394  C* SUBROUTINE USED FOR FIELD LINE TRACIN Line 394  C* SUBROUTINE USED FOR FIELD LINE TRACIN
394  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
395  C*******************************************************************  C*******************************************************************
396        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
397        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
398  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
399        ZM=P(3)                                                                  ZM=P(3)                                                          
400        FLI=P(1)*P(1)+P(2)*P(2)+1E-15        FLI=P(1)*P(1)+P(2)*P(2)+1E-15
# Line 476  C                 TO THE LOCAL GEODETIC Line 476  C                 TO THE LOCAL GEODETIC
476  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
477  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
478  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
479        DIMENSION         V(3),B(3),G(144)        DIMENSION         V(3),B(3),G(196)
480        CHARACTER*258      NAME        CHARACTER*258      NAME
481        INTEGER NMAX        INTEGER NMAX
482        REAL TIME        REAL TIME
483        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
484        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
485        SAVE/MODEL/        SAVE/MODEL/
486        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 581  C--------------------------------------- Line 581  C---------------------------------------
581        CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
582        CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
583  C     ### FILMOD, DTEMOD arrays +1  C     ### FILMOD, DTEMOD arrays +1
584  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
585        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)
586        DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
587        INTEGER L1,L2,L3        INTEGER L1,L2,L3
588        INTEGER NMAX        INTEGER NMAX
# Line 607  c     COEF3 = COEFPATH(1:16)//COEF3 Line 607  c     COEF3 = COEFPATH(1:16)//COEF3
607  c     FILMOD(1) = COEF1  c     FILMOD(1) = COEF1
608  c     FILMOD(2) = COEF2  c     FILMOD(2) = COEF2
609  c     FILMOD(3) = COEF3  c     FILMOD(3) = COEF3
610  c      print *, "qui"        print *, "qui"
611        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
612        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
613        FILMOD(3) = P3(1:L3)        FILMOD(3) = P3(1:L3)
614  c      print *, "qua"        print *, "qua"
615  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
616  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
617  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
618  c     WRITE(*,*) FILMOD(1)        WRITE(*,*) FILMOD(1)
619  c     WRITE(*,*) FILMOD(2)        WRITE(*,*) FILMOD(2)
620  c     WRITE(*,*) FILMOD(3)        WRITE(*,*) FILMOD(3)
621  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
622        DATA   DTEMOD / 2000., 2005., 2010./        DATA   DTEMOD / 2005., 2010., 2015./
623  c      c    
624  c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',            
625  c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',      
# Line 633  C     ### numye = numye + 1 ; is number Line 633  C     ### numye = numye + 1 ; is number
633  C      C    
634  c     NUMYE=13  c     NUMYE=13
635        NUMYE=2        NUMYE=2
636  c      print *, "quo"        print *, "quo"
637                
638  C      C    
639  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
# Line 652  c     L = (IYEA - 1945)/5 + 1 Line 652  c     L = (IYEA - 1945)/5 + 1
652        FIL1 = FILMOD(L)          FIL1 = FILMOD(L)  
653        DTE2 = DTEMOD(L+1)        DTE2 = DTEMOD(L+1)
654        FIL2 = FILMOD(L+1)        FIL2 = FILMOD(L+1)
655  c      print *, "que"        WRITE(*,*) FIL1
656          WRITE(*,*) FIL2
657          print *, "que"
658  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
659        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
660        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
661          print *, "quessss"
662        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
663        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
664  c      print *, "quj"        print *, "quj"
665  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
666        IF (L .LE. NUMYE-1) THEN                                IF (L .LE. NUMYE-1) THEN                        
667           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
# Line 667  C--   DETERMINE IGRF COEFFICIENTS FOR YE Line 670  C--   DETERMINE IGRF COEFFICIENTS FOR YE
670           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,               CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
671       1        GH2, NMAX, GHA)                                           1        GH2, NMAX, GHA)                                    
672        ENDIF        ENDIF
673  c      print *, "quw"        print *, "quw"
674  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
675        F0=0.D0        F0=0.D0
676        DO 1234 J=1,3        DO 1234 J=1,3
# Line 741  C     ---------------------------------- Line 744  C     ----------------------------------
744        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
745  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
746   667  FORMAT(A258)   667  FORMAT(A258)
747  c      print *," gui"        print *," gui"
748        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
749  c      print *," gua"        print *," gua"
750        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
751  c      print *," gue"        print *," gue"
752        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
753  c      print *," guo"        print *," guo"
754  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
755  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
756  C                                                                                C                                                                              
# Line 784  C -------------------------------------- Line 787  C --------------------------------------
787                  ENDIF                                                                  ENDIF                                                
788  2233        CONTINUE                                                      2233        CONTINUE                                                    
789  2211    CONTINUE                                                          2211    CONTINUE                                                        
790  c      print *," guj"        print *," guj"
791                                                                                                                                                                    
792  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
793  c      print *," guw IER",IER        print *," guw IER",IER
794        if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions        if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions
795                                                                                                                                                                    
796          RETURN                                                                RETURN                                                      
# Line 944  C ERA, EREQU and ERPOL as recommended by Line 947  C ERA, EREQU and ERPOL as recommended by
947  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
948  C-----------------------------------------------------------------  C-----------------------------------------------------------------
949          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
950          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=258) TP1,TP2,TP3
951          INTEGER L1,L2,L3          INTEGER L1,L2,L3
952          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
953          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23