/[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.6 by mocchiut, Tue May 15 14:31:14 2012 UTC revision 1.12 by mocchiut, Mon Jan 19 12:32:13 2015 UTC
# Line 1  Line 1 
         subroutine igrf_sub(xlat,xlong,year,height,  
      &          xl,icode,dip,dec)  
 c----------------------------------------------------------------  
 c   INPUT:  
 c     xlat     geodatic latitude in degrees  
 c     xlong    geodatic longitude in degrees  
 c     year     decimal year (year+month/12.0-0.5 or year+day-of-year/365  
 c          or 366 if leap year)  
 c     height     height in km  
 c   OUTPUT:  
 c     xl     L value  
 c     icode     =1  L is correct; =2  L is not correct;  
 c          =3  an approximation is used  
 c     dip     geomagnetic inclination in degrees  
 c     dec     geomagnetic declination in degress  
 c----------------------------------------------------------------  
   
       REAL              LATI,LONGI  
       COMMON/GENER/     UMR,ERA,AQUAD,BQUAD  
       SAVE /GENER/  
1  C  C
       CALL INITIZE  
       ibbb=0  
       ALOG2=ALOG(2.)  
       ISTART=1  
         lati=xlat  
         longi=xlong  
 c  
 C----------------CALCULATE PROFILES-----------------------------------  
 c  
         CALL FELDCOF(YEAR,DIMO)  
         CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)  
         CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)  
         DIP=ASIN(BDOWN/BABS)/UMR  
       DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR  
       RETURN  
       END  
 c  
 c  
2  C SHELLIG.FOR, Version 2.0, January 1992  C SHELLIG.FOR, Version 2.0, January 1992
3  C  C
4  C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2    C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2  
# Line 196  C--------------------------------------- Line 158  C---------------------------------------
158        SAVE /FIDB0/        SAVE /FIDB0/
159        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
160        SAVE /GENER/        SAVE /GENER/
161          REAL FLS
162  C  C
163  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
164  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING
# Line 203  C-- STEQ IS STEP SIZE FOR INTEGRATION Line 166  C-- STEQ IS STEP SIZE FOR INTEGRATION
166  C  C
167        DATA RMIN,RMAX    /0.05,1.01/        DATA RMIN,RMAX    /0.05,1.01/
168        DATA STEP,STEQ    /0.20,0.03/        DATA STEP,STEQ    /0.20,0.03/
169          BEQU=1.E10        BEQU=1.E10
170          FLS = FL
171  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES
172        RLAT=GLAT*UMR        RLAT=GLAT*UMR
173        CT=SIN(RLAT)                                                      CT=SIN(RLAT)                                              
# Line 304  C*****INNER LOOP (FOR QUADRATURE)       Line 268  C*****INNER LOOP (FOR QUADRATURE)      
268        HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                                    HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                            
269        ZQ=Z*Z        ZQ=Z*Z
270        R=HLI+SQRT(HLI*HLI+ZQ)        R=HLI+SQRT(HLI*HLI+ZQ)
271          IF(R.NE.R)THEN
272             FL = FLS
273             RETURN
274          ENDIF
275        IF(R.LE.RMIN)GOTO30                                      IF(R.LE.RMIN)GOTO30                              
276        RQ=R*R        RQ=R*R
277        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
# Line 580  C     ### updated to IGRF-2005 version - Line 548  C     ### updated to IGRF-2005 version -
548  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
549        CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
550        CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
551  C     ### FILMOD, DTEMOD arrays +1        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(2)
 c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  
       DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)  
552        DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
553        INTEGER L1,L2,L3        INTEGER L1,L2,I1
554        INTEGER NMAX        INTEGER NMAX
555        REAL TIME        REAL TIME
556        CHARACTER *258 P1,P2,P3        CHARACTER *258 P1,P2
557        COMMON/PPATH/ L1,L2,L3,P1, P2, P3        COMMON/PPATH/ I1,L1,L2,P1,P2
558        SAVE/PPATH/        SAVE/PPATH/
559        COMMON/MODEL/   GH1,NMAX,TIME,FIL1        COMMON/MODEL/   GH1,NMAX,TIME,FIL1
560        SAVE/MODEL/        SAVE/MODEL/
561        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD
562        SAVE/GENER/        SAVE/GENER/
563  C     ### updated to 2005  c      print *, "qui"
 C     CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80  
         
 c     COEFPATH  = 'OrbitalInfo/src/'  
 c     COEF1 = 'dgrf00.dat'  
 c     COEF2 = 'igrf05.dat'  
 c     COEF3 = 'igrf05s.dat'  
 c     COEF1 = COEFPATH(1:16)//COEF1  
 c     COEF2 = COEFPATH(1:16)//COEF2  
 c     COEF3 = COEFPATH(1:16)//COEF3  
 c     FILMOD(1) = COEF1  
 c     FILMOD(2) = COEF2  
 c     FILMOD(3) = COEF3  
       print *, "qui"  
564        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
565        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
566        FILMOD(3) = P3(1:L3)  c      print *, "qua"
       print *, "qua"  
 c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  
 c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  
 c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  
       WRITE(*,*) FILMOD(1)  
       WRITE(*,*) FILMOD(2)  
       WRITE(*,*) FILMOD(3)  
 c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  
       DATA   DTEMOD / 2005., 2010., 2015./  
 c      
 c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              
 c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        
 c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',        
 c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',  
 c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/  
 c     DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,  
 c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./        
 C      
 C     ### numye = numye + 1 ; is number of years represented by IGRF  
 C      
 c     NUMYE=13  
       NUMYE=2  
       print *, "quo"  
         
567  C      C    
568  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
569  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
# Line 644  C     Line 573  C    
573  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR
574        TIME = YEAR        TIME = YEAR
575        IYEA = INT(YEAR/5.)*5        IYEA = INT(YEAR/5.)*5
576  c     L = (IYEA - 1945)/5 + 1        L = IYEA + 5
577        L = (IYEA - 2000)/5 + 1  C
578        IF(L.LT.1) L=1        DTE1 = REAL(IYEA)  
579        IF(L.GT.NUMYE) L=NUMYE                FIL1 = FILMOD(1)  
580        DTE1 = DTEMOD(L)          DTE2 = REAL(L)
581        FIL1 = FILMOD(L)          FIL2 = FILMOD(2)
582        DTE2 = DTEMOD(L+1)  c      print *,'IYEA ',IYEA,' L ',L
583        FIL2 = FILMOD(L+1)  c      WRITE(*,*) FIL1
584        WRITE(*,*) FIL1  c      WRITE(*,*) FIL2
585        WRITE(*,*) FIL2  c      print *, "que"
       print *, "que"  
586  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
587        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
588        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
589        print *, "quessss"  c      print *, "quessss"
590        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
591        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
592        print *, "quj"  c      print *, "quj"
593  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
594        IF (L .LE. NUMYE-1) THEN                                IF (I1.EQ.0) THEN      
595           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
596       1        NMAX2, GH2, NMAX, GHA)                               1        NMAX2, GH2, NMAX, GHA)                        
597        ELSE                      ELSE              
598           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,               CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
599       1        GH2, NMAX, GHA)                                           1        GH2, NMAX, GHA)                                    
600        ENDIF        ENDIF
601        print *, "quw"  c      print *, "quw"
602  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
603        F0=0.D0        F0=0.D0
604        DO 1234 J=1,3        DO 1234 J=1,3
605           F = GHA(J) * 1.D-5           F = GHA(J) * 1.D-5
606           F0 = F0 + F * F           F0 = F0 + F * F
607   1234 CONTINUE   1234 CONTINUE
608        DIMO = DSQRT(F0)        DIMO = REAL(DSQRT(F0))
609                
610        GH1(1) =  0.0        GH1(1) =  0.0
611        I=2                  I=2          
# Line 693  c      print *, "quq" Line 621  c      print *, "quq"
621           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
622           F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
623           IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
624           GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
625           I = I+1                                                   I = I+1                                        
626           DO 9 M=1,N                                               DO 9 M=1,N                                    
627              F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
628              IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))              IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))
629              GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
630              GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
631              I=I+2              I=I+2
632   9       CONTINUE                                             9       CONTINUE                                          
633        RETURN        RETURN
# Line 744  C     ---------------------------------- Line 672  C     ----------------------------------
672        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
673  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
674   667  FORMAT(A258)   667  FORMAT(A258)
675        print *," gui"  c      print *," gui"
676        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
677        print *," gua"  c      print *," gua"
678        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
679        print *," gue"  c      print *," gue"
680        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
681        print *," guo"  c      print *," guo"
682  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
683  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
684  C                                                                                C                                                                              
# Line 787  C -------------------------------------- Line 715  C --------------------------------------
715                  ENDIF                                                                  ENDIF                                                
716  2233        CONTINUE                                                      2233        CONTINUE                                                    
717  2211    CONTINUE                                                          2211    CONTINUE                                                        
718        print *," guj"  c      print *," guj"
719                                                                                                                                                                    
720  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
721        print *," guw IER",IER  c      print *," guw IER",IER
722        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
723                                                                                                                                                                    
724          RETURN                                                                RETURN                                                      
# Line 931  C -------------------------------------- Line 859  C --------------------------------------
859          END                                                                      END                                                            
860  C  C
861  C  C
862          SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3)          SUBROUTINE INITIZE(ISSEC,TP1,TL1,TP2,TL2)
863  C----------------------------------------------------------------  C----------------------------------------------------------------
864  C Initializes the parameters in COMMON/GENER/  C Initializes the parameters in COMMON/GENER/
865  C  C
# Line 946  C Line 874  C
874  C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL  C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL
875  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
876  C-----------------------------------------------------------------  C-----------------------------------------------------------------
877          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,ISSEC
878          CHARACTER (len=258) TP1,TP2,TP3          CHARACTER (len=*) :: TP1,TP2
879          INTEGER L1,L2,L3          INTEGER L1,L2
880          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2
881          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ I1,L1,L2,P1,P2
882          SAVE/PPATH/          SAVE/PPATH/
883    
884          COMMON/GENER/UMR,ERA,AQUAD,BQUAD          COMMON/GENER/UMR,ERA,AQUAD,BQUAD
885          SAVE/GENER/          SAVE/GENER/
886    
887            I1 = ISSEC
888          L1 = TL1          L1 = TL1
889          L2 = TL2          L2 = TL2
890          L3 = TL3  
891          P1 = TP1(1:L1)          P1 = TP1(1:L1)
892          P2 = TP2(1:L2)          P2 = TP2(1:L2)
         P3 = TP3(1:L3)  
893    
894          ERA=6371.2          ERA=6371.2
895          EREQU=6378.16          EREQU=6378.16
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23