--- /dev/null
+C Last update: 8-4-1993 ore 17:20
+C - Version 0.0 - Debbugging in progress
+C Stato attuale:
+C - Questa versione ha un disco e 2 sfere o roche,
+C senza grafica, limb darkening non testato, senza common envelope.
+C Debbugging in corso.
+C -------------------------------------------------------------
+C History:
+C february 1993 : performing test 8 - test 9
+C 1-3-1993 : changed omega definition to omega-q**2/2/(1+q)
+C -------------------------------------------------------------
+C ***********************************************************
+C
+C CCCCCC BBBBBB SSSSS
+C CC BB B S
+C CC BBBBBB SS
+C CC BB B S
+C CCCCCC BBBBBB SSSSS
+C
+C
+C C L O S E B I N A R Y S Y S T E M
+C
+C L I G T H C U R V E A N A L Y S I S P R O G R A M
+C
+C ***********************************************************
+C Marcello Galli - Loretta Solmi
+C ***********************************************************
+C
+C Version 0.0 : 14-5-1991 : writing in progress
+C
+C ***********************************************************
+C
+C Given the parameters of the bynary system, this program
+C computes the light curve for given colors, accounting
+C for limb darkening, gravity darkening and reflection
+C effects.
+C It is possible to simulate the presence of a disk around
+C one component of the system.
+C
+C *****************************************************************
+C Units: CGS,
+C For Roche model the unit distance is the distance between stars A-B
+C Intensity : erg/sec/cm**2/sterad ,
+C F = ca/4pi T**4 * area = erg/sec/sterad
+C with c*a/(4*pi)=1.8065 erg/sec/cm**2/K**4/sterad
+C to shorten calculations i should use adimensional units.
+C -----------------------------------------------------------------
+C Note per la messa a punto:
+C - Ci sono ridondanze nel common /allin/:
+C Rxij,Ryij,Rzij, iaus usati solo fino a call CALCALL escluso;
+C Xcosgi,fsour usati solo da CALCALL compresa in poi.
+C - Qualcosa in LUCE1 puo' essere messo integer*2 (IGAINED,ilost)
+C - Il problema principale sono i tempi di esecuzione della riflessione
+C sto cercando di usare un doppio grid (fine per eclisse, coarse per
+C riflessione) in modo da abbreviare il tempo di calcolo della
+C riflessione. I path sono moltissimi, per una stella descritta
+C da un 10000 surf.elem.,ove eclisse e' precisa piu' di qualche %,
+C si va sul mezzo milione di paths,ci si sta a fatica anche come spazio.
+C Poiche' la riflessione per 2 stelle eguali a contatto va sul
+C 3-4 % della luce dell'oggetto, diminuire il numero di punti
+C per la riflessione di un ordine di grandezza potrebbe andare.
+C Le subroutine di riflessione hanno tempo ed ingombro accettabile
+C per oggetti descritti da circa un 100 surface elements.
+C Altro modo di ridurre i punti e' usare surf. elements con
+C egual area uinvece che egual numero di fi per ogni theta.
+C Questo riduce il numero di punti di circa 1/3, a parita'
+C di precisione.
+C - L'inserimento in un metodo di minimo totale puo' consistere
+C nel trasformare il main , dal 300 continue in poi,
+C in una routine, chiamata dalla subroutine di minimo, che
+C lo pilota tramite i flags e parametri, a seconda di cosa cambia
+C nel modello; ma i parametri e flags non sono ben pensati
+C per questo uso.
+C - Benche' parte del programma sia predisposta per un calcolo
+C su diverse bande di luce, la riflessione tratta per ora
+C solo il bolometrico.
+C qui non mi e' chiaro come procedere; ora con la riflessione
+C potrei (ma non ancora c'e') aumentare o T e poi ricalcolare
+C i flussi della banda; oppure,
+C nel caso in cui do' in input il flusso di ogni banda,
+C divido il flusso bol.+riflesso bol. nelle bande secondo l'input.
+C PER QUANTO RIGUARDA L'ALBEDO abbiamo una ALBEDO eguale per
+C ogni oggetto (PAR(ippar(9))) che e' semplicemente un flag, che
+C andando da 0 ad 1 attenua l'effetto della riflessione ed
+C una albedo (PARS(12,IS)), definita per ogni oggetto,
+C che fa la stessa cosa. Potra' essere comodo in input, ma
+C averne 2 intruduce un conto in piu' in un punto critico del programma.
+C Conviene fare una riflessione con albedo dipendente dall'oggetto
+C e dal colore? Questo equivale ad un modello
+C di atmosfera e non ne so molto.
+C - Gli elementi di superficie sono in una unica lista,
+C indipendentemente dalla stella cui appartengono, a parte
+C il risparmio di spazio, serve a qualcosa cio', visto che
+C spesso tratto oggetto per oggetto e devo portarmi dietro
+C una collezione di puntatori per sapere dove iniziano e finiscono?
+C - Spazi ausiliari probabilmente ridondanti:IAUS,NUMOBJ,AUS,AUS1
+C - GC e' usato solo per rinormalizzare GCX,GCY,GCZ, si puo' non conservare
+C - Bisogna inserire controlli sulla possibilita' di avere
+C oggetti intersecantisi. Il programma non se ne accorge e pasticcia
+C riflessione ed eclissi.
+C - Approssimazioni numeriche nei sin e cos fanno si che aree e
+C distanze abbiano solo 4 cifre significative, qualcosa forse
+C andra' messo in doppia precisione. In certi punti si rimedia a
+C questo mettendo =0 le cose sotto R4prec. Ma, per esempio in CYL,
+C thin disk a spessore varabile, quando si usano sui 50000 punti
+C le approssimazioni numeriche danno gia' nel calcolo area tot
+C un errore circa del 0,04 %
+C - Per come e' ora ALLSEL dimcx,y,z puo' essere abs(dimcx,y,z)
+C con ovvio risparmio di tempo.
+C - Le routines lobesf1,lobesf2,lobesf2 sono da modificare,
+C riscrivendo le espressioni in modo da eliminare le potenze
+C - usare adimensional units per accorciare calcolo flussi da T**4
+C mettendo: ac/4pi=AC4PI=1
+C - ridondanze ed ineleganze in opzioni di stampe ed altre.
+C la frazione di luce in una banda per un oggetto in PAR
+C andrebbe usata per metterci il valore calcolato , ora in
+C /filtri/ .. compfrac(banda,object). Non si puo' finche'
+C e' usata anche per segnalare che T usare nella plankiana.
+C Anche i vari modflag delle bande di colore possono non essere OK.
+C - Bisogna notare che il coefficiente di gravity darkening e',
+C seguendo Wilson, g=4*beta. Beta e' la notazione di Lucey
+C - Per finire: l'inglese dei commenti e' lievemente penoso.
+C
+C ***********************************************************
+C ARRAYS AND VARIABLES DESCRIPTION
+C
+C PARAMETER Statements are used for compile-time dimensioning:
+C MAXPT= number of fine grid surface points describing the 3 objects
+C MAXPTC= number of coarse grid surface points describing the 3 objects
+C MAXALLIN= number of reflection paths
+C MAXFASI= number of values for the phase
+C maxbande= 5 usually UBVRI
+C maxflags= number of flags defining the program run
+C maxfval=number of flag value descriptors
+C maxpar= number of parameters describing the system
+C maxpars= number of parameters describing each object
+C MAXCELL= number of grid points used in light integration
+C maxtitle= number of max title line for printed output
+C maxfilt= number of max colour bands
+C maxlam= number of max lambda points describing a colour band
+C
+C COMMON/SURF/ fine grid description the surface elements of the system:
+C NPNTMX=MAXPT=maximum number of surface elements used to describe
+C the 3 objects of the system: the two stars and the disk.
+C NPNT=actual number of surface elements
+C NPNT1= number of surface elements for object 1
+C NPNT2,NPNT3= same for objects 2 and 3
+C N1I,N1F= first and last surface element for object 1 in following arrays
+C N2I,N2F,N3I,N3F= same for objects 2 and 3
+C X,Y,Z= position of each surface element in the corotating system
+C G= surface gravity
+C GX,GY,GZ= component of the normal unit vector to the surface element
+C FX,FY,FZ= components of the surf.element tangent vector (along phi)
+C TX,TY,TZ= components of the surf.element tangent vector (along theta)
+C T= temperature of the surface element
+C A= area of the surface element
+C FB=bolometric flux of the surface element (reflection excluded)
+C F= flux for each color
+C FRIFL= flux reflected from the surface element
+C NUMOBJ= number of the object for each surface element
+C XLIMB= limb darkening coefficient for each surf. element
+C ALB= albedo for each surf. element
+C
+C COMMON/SURFC/ describes the surface elements of the system:
+C similar to /SURF/, but coarse grid, used for reflection computation.
+C DIMCX,DIMCY,DIMCZ : lenght of the coarse surf.element
+C
+C COMMON/ALLIN/ contains data used to compute reflection effects,
+C namely a list of reflection paths (alignement's list) containig:
+C IAUS is used as a logical existence flag in subroutine ALLSEL
+C
+C COMMON/LIGHT/ contains phases, observed and computed light
+C for each phase and color. maxbande of this common is the
+C same as maxfilt of common filtri, this duplication allow
+C a more easy management of common dimension during test phases
+C FASE= phase values
+C O = observed flux for each colour band and phase
+C C= computed flux for each band and object and phase
+C CT= C summed on objects
+C CBOL = bolometric flux for each object and phase
+C CBOLT = CBOL summen ob objects
+C CBOLR = bolometric reflected flux for each object and phase
+C CBOLT = CBOLR summen ob objects
+C AREA= total surface area of the star
+C FBOLTOT = flux norm factor : ( L norm.*FBOLTOT = true L )
+C FBOL= bolometric tot. luminosity flux for the object without reflected
+C FBOLREF= bolometric tot luminosity reflected by the object
+C
+C COMMON/VISUAL/ contains an array used to compute the light
+C received by the observed, namely a grid of points on the
+C plane perpendicular to the observer over which the image
+C of the system is proiected
+C
+C COMMON /ROCHE/ roche model Lagrange points and lobe'e dimension
+C AL1,AL2,AL3 = x coord of lagrange points : 1,2,3
+C XA1,XA2 = intersection of star A surface with X axis
+C XB1,XB2 = same for star B (XB1,XA1 are near inner Lagr.point)
+C XB2 near L2, XA2 near L3
+C YA,YB = approx transverse dimesion of critical lobes(common L1)
+C See: Plavec - BAC - Vol 15 :165 (1964)
+C omegal1,omegal2,omegal3 = omega for critical lobes
+C rpole,xpole,ypole,zpole,tpole,gxpole,gypole,gzpole = position, radius,
+C temperature, surf. versor for pole surf.element
+C
+C COMMON/FILTRI/ contains the colour bands description, filled in
+C routine FILFILT. maxfilt of this common is the
+C same as maxbande of common light , this duplication allow
+C a more easy management of common dimension during test phases
+C NLAMBDA = number of lambda value describing the band
+C deltalam = constant step for lambda values
+C alambda= lambda values (Anstrong)
+C trasmiss= filter trasmission coefficients for the corresponding lambda
+C compfrac=computed fraction of intensity for each object in each band
+C
+C COMMON/PAR/ contains the calculation parameters, which
+C can be modified by the user input
+C
+C COMMON/TITLE/ some title lines to be printed on output
+C
+C -----------------------------------------------------------------
+C
+ PARAMETER (MAXPT=200000 )
+ PARAMETER (MAXPTC=50000 ,MAXALLIN=250000)
+ PARAMETER (MAXFASI=1000, MAXBANDE=5)
+ PARAMETER (MAXFLAG=21 , MAXPAR=128 , MAXPARS=20 )
+ PARAMETER (MAXCELL=1000 )
+ PARAMETER (MAXTITLE=5)
+ PARAMETER ( MAXFILT=5 , MAXLAM=15 )
+C
+ COMMON /SURF/ NPNTMX,NPNT,NPNT1,NPNT2,NPNT3,
+ A N1I,N1F,N2I,N2F,N3I,N3F,
+ 1 X(MAXPT),Y(MAXPT),Z(MAXPT),
+ 2 G(MAXPT),GX(MAXPT),GY(MAXPT),GZ(MAXPT),
+ 3 FX(MAXPT),FY(MAXPT),FZ(MAXPT),
+ 4 TX(MAXPT),TY(MAXPT),TZ(MAXPT),
+ 5 T(MAXPT),A(MAXPT),
+ 6 FB(MAXPT),F(MAXPT,MAXBANDE),FRIFL(MAXPT),
+ 7 XLIMB(MAXPT),NUMOBJ(MAXPT),ALB(MAXPT),
+ 8 NUMCOARSE(MAXPT)
+ DIMENSION NP123(3)
+ DIMENSION NIF(2,3)
+ EQUIVALENCE (NP123(1),NPNT1),(NIF(1,1),N1I)
+C
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+ COMMON /ALLIN/ NALLMX,NALL,ITERDONE,
+ 1 ISOUR(MAXALLIN),JRIC(MAXALLIN),
+ 2 TRANSFI(MAXALLIN),TRANSFJ(MAXALLIN),
+ 3 FSOURI(MAXALLIN),FSOURJ(MAXALLIN),
+ 4 FRICI(MAXALLIN),FRICJ(MAXALLIN),
+ 5 RIJ(MAXALLIN),
+ 6 RXIJ(MAXALLIN),RYIJ(MAXALLIN),RZIJ(MAXALLIN),
+ 7 COSGI(MAXALLIN),COSGJ(MAXALLIN),
+ 8 IAUS(MAXALLIN)
+C
+ COMMON /LIGHT/ NFASIMX,NFASI,NBANDE,
+ 1 FASE(MAXFASI),O(MAXBANDE,MAXFASI),
+ 2 C(3,MAXBANDE,MAXFASI),CT(MAXBANDE,MAXFASI),
+ 3 CBOL(3,MAXFASI),CBOLT(MAXFASI),
+ 4 CBOLR(3,MAXFASI),CBOLTR(MAXFASI),
+ 5 AREA(3),FBOLTOT(3),FBOL(3),FBOLREF(3)
+C
+ COMMON /VISUAL/ NCELL,CELLMAX,CELLMIN,DCELL,
+ 1 XCELL(MAXCELL),YCELL(MAXCELL)
+ DIMENSION ICELL(MAXCELL,MAXCELL)
+C
+ COMMON /ROCHE/ AL1,AL2,AL3,XA1,XA2,XB1,XB2,YA,YB,
+ 1 OMEGAL1,OMEGAL2,OMEGAL3,
+ 2 RPOLE(3),XPOLE(3),YPOLE(3),ZPOLE(3),TPOLE(3),
+ 3 GXPOLE(3),GYPOLE(3),GZPOLE(3)
+C
+ COMMON /FILTRI/ NFILT,NLAM, NLAMBDA(MAXFILT),DELTALAM(MAXLAM),
+ 1 ALAMBDA(MAXLAM,MAXFILT),TRASMISS(MAXLAM,MAXFILT),
+ 2 COMPFRAC(MAXFILT,3)
+C
+ COMMON /PAR/ NFLAG,IFLAG(MAXFLAG),NPAR,PAR(MAXPAR)
+ DIMENSION PARS(MAXPARS,3)
+ EQUIVALENCE (PARS(1,1),PAR(1))
+C
+ COMMON/TITLE/NTITLEMX,NTITLE,TITLE(MAXTITLE)
+ CHARACTER*80 TITLE
+C
+ DIMENSION IPPAR1(MAXFLAG),IPPAR2(MAXFLAG),MODFLAG(MAXFLAG)
+C
+ CHARACTER *10 NAMFLAG(MAXFLAG)
+ CHARACTER *20 DESCFLAG(MAXFLAG)
+C
+ CHARACTER *10 NAMPAR(MAXPAR)
+ CHARACTER *20 DESCPAR(MAXPAR)
+C
+ DIMENSION IFLAGDEF(MAXFLAG),PARDEF(MAXPAR)
+C
+ DIMENSION AUS(MAXPT),AUS1(MAXPT),AUS2(MAXPT),AUS3(MAXPT)
+C aus,aus1 used by subroutine sphere (where dimensioned nfi)
+C aus,aus1,aus2,aus3 are used by subroutine luce
+C (where dimensioned maxpt)
+C aus used in subroutine CYL (where dimensioned NTHF*5; NTHF=coarse factor)
+C aus used in subroutine PRINTING (dim(MAXPT) to store radii to print
+C
+C -------------------------------------------------------------
+ DATA PI /3.1415926/
+C a*c/4/pi , flux=ac4pi T**4 = erg/cm**2/sec/sterad
+ DATA AC4PI /1.8065E-5/
+C
+ LOGICAL KTEST,STOP
+ LOGICAL PRINTFILE,PRINT6
+ LOGICAL PLOTFILE,PLOT6
+ EXTERNAL ROCHESF1,ROCHESF2
+C
+C The following data define defaults for parameters and flags
+C along with a name and a descriptor for each of them.
+C For each flag ippar1 and ippar2 are two pointers to the
+C parameters, which identify the range of parameters belonging
+C to the flag. Same parameters are hidden i.e. not referred to.
+C
+C --------------------------------------------------------------
+ DATA (NAMFLAG(I),DESCFLAG(I),IFLAGDEF(I),
+ A IPPAR1(I),IPPAR2(I),I=1,MAXFLAG)
+ 1 /'STARA ',' Star A ', 1, 1,20,
+ 2 'STARB ',' Star B ', 0, 21,40,
+ 3 'DISK ',' Disk ', 0, 41,60,
+ 4 'U ',' U-color band ', 0, 61,65,
+ 5 'B ',' B-color band ', 0, 69,73,
+ 6 'V ',' V-color band ', 0, 77,81,
+ 7 'R ',' R-color band ', 0, 85,89,
+ 8 'I ',' I-color band ', 0, 93,97,
+ 9 'REFLECTION',' Reflection comput. ', 1, 101,104,
+ A 'ZEROREFL ',' Zero refl. source ', 0, 0, 0,
+ 1 'ORBIT ',' Orbital parameters ', 0, 105,108,
+ 2 'PHASES ',' Sets phases ', 1, 111,113,
+ 3 'GRID ',' Sets project.grid', 1, 114,116,
+ 4 'PARAMETER ',' Sets a parameter ', 0, 1, MAXPAR,
+ 5 'GO ',' Program runs ', 0, 0, 0,
+ 6 'STOP ',' Program stops ', 0, 0, 0,
+ 7 'EXIT ',' Program stops ', 0, 0, 0,
+ 8 'FLAGS ',' Sets a flag ', 0, 0, 0,
+ 9 'OFF ',' Set next flags off ', 0, 0, 0,
+ A 'PRINT ',' Output is printed ', 1, 117,122,
+ B 'PLOT ',' Plot produced ', 0, 123,127 /
+C
+ DATA (NAMPAR(I),DESCPAR(I),PARDEF(I),I=1,20)
+ 1 /'SHAPEA ' , ' sphere=1,roche=2 ' , 2.0,
+ 2 'RA ' , ' radius of star A ' , 0.0,
+ 3 'X0A ' , ' X position of A ' , 0.0,
+ 4 'Y0A ' , ' Y position of A ' , 0.0,
+ 5 'Z0A ' , ' Z position of A ' , 0.0,
+ 6 'OMEGAA ' , ' potential of A ' , 0.0,
+ 7 'MESHA ' , ' num. of theta mesh ' , 5.0,
+ 8 'NPHIA ' , ' unused ' , 0.0,
+ 9 'GA ' , ' gravity dark. of A ' , 0.0,
+ O 'TEMPA ' , ' pole temperature A ' , 1.0,
+ 1 'BOLA ' , ' bolometric lum.of A' , 0.0,
+ 2 'ALBA ' , ' bolometric albedo A' , 1.0,
+ 3 'LIMBA ' , ' limb darkening of A' , 0.0,
+ 4 'CORDA ' , ' arch,cord,tang.seg.' , 3.0,
+ 5 'RIA ' , ' disk inner radius ' , 0.0,
+ 6 'HA ' , ' disk height at R ' , 0.0,
+ 7 'HIA ' , ' disk height at RI ' , 0.0,
+ 8 'INCLINA ' , ' ang.z-z orbit.plane' , 0.0,
+ 9 'ROTATEA ' , ' ang.x-x orbit plane' , 0.0,
+ O 'INNERA ' , ' unused ' , 0.0 /
+ DATA (NAMPAR(I),DESCPAR(I),PARDEF(I),I=21,40)
+ 1 /'SHAPEB ' , ' sphere=1,roche=2 ' , 1.0,
+ 2 'RB ' , ' radius of star B ' , 0.0,
+ 3 'X0B ' , ' X position of B ' , 0.0,
+ 4 'Y0B ' , ' Y position of B ' , 0.0,
+ 5 'Z0B ' , ' Z position of B ' , 0.0,
+ 6 'OMEGAB ' , ' potential of B ' , 0.0,
+ 7 'MESHB ' , ' num. of theta mesh ' , 5.0,
+ 8 'NPHIB ' , ' unused ' , 0.0,
+ 9 'GB ' , ' gravity dark. of B ' , 0.0,
+ O 'TEMPB ' , ' pole temperature B ' , 1.0,
+ 1 'BOLB ' , ' bolometric lum.of B' , 0.0,
+ 2 'ALBB ' , ' bolometric albedo B' , 1.0,
+ 3 'LIMBB ' , ' limb darkening of B' , 0.0,
+ 4 'CORDB ' , ' arch,cord,tang.seg.' , 3.0,
+ 5 'RIB ' , ' disk inner radius ' , 0.0,
+ 6 'HB ' , ' disk height at R ' , 0.0,
+ 7 'HIB ' , ' disk height at RI ' , 0.0,
+ 8 'INCLINB ' , ' ang.z-z orbit.plane' , 0.0,
+ 9 'ROTATEB ' , ' ang.x-x orbit plane' , 0.0,
+ O 'INNERB ' , ' unused ' , 0.0 /
+ DATA (NAMPAR(I),DESCPAR(I),PARDEF(I),I=41,60)
+ 1 /'SHAPEC ' , ' disk=3 ' , 3.0,
+ 2 'RC ' , ' radius of disk:C ' , 4.0,
+ 3 'X0C ' , ' X position of disk ' , 0.0,
+ 4 'Y0C ' , ' Y position of disk ' , 0.0,
+ 5 'Z0C ' , ' Z position of disk ' , 0.0,
+ 6 'OMEGAC ' , ' potential of C ' , 0.0,
+ 7 'MESHC ' , ' num. of theta mesh ' , 4.0,
+ 8 'NPHIC ' , ' unused ' , 0.0,
+ 9 'GC ' , ' gravity dark. of C ' , 0.0,
+ O 'TEMPC ' , ' pole temperature C ' , 1.0,
+ 1 'BOLC ' , ' bolometric lum.of C' , 0.0,
+ 2 'ALBC ' , ' bolometric albedo C' , 1.0,
+ 3 'LIMBC ' , ' limb darkening of C' , 0.0,
+ 4 'CORDC ' , ' arch,cord,tang.seg.' , 0.0,
+ 5 'RIC ' , ' disk inner radius ' , 0.0,
+ 6 'HC ' , ' disk height at R ' , 0.0,
+ 7 'HIC ' , ' disk height at RI ' , 0.0,
+ 8 'INCLINC ' , ' ang.z-z orbit.plane' , 0.0,
+ 9 'ROTATEC ' , ' ang.x-x orbit plane' , 0.0,
+ O 'INNERC ' , ' inner surf drawn ' , 0.0 /
+C
+ DATA (NAMPAR(I),DESCPAR(I),PARDEF(I),I=61,100)
+ 1 /'ULAMB1 ' , ' U color lambda 1 ' , 3300.0,
+ 2 'ULAMB2 ' , ' U color lambda 2 ' , 3700.0,
+ 3 'ULUM_A ' , ' frac. U lumin.for A' , 0.0,
+ 4 'ULUM_B ' , ' frac. U lumin.for B' , 0.0,
+ 5 'ULUM_C ' , ' frac. U lumin.for C' , 0.0,
+ 6 'UALB_A ' , ' albedo U for A ' , 1.0,
+ 7 'UALB_B ' , ' albedo U for B ' , 1.0,
+ 8 'UALB_C ' , ' albedo U for C ' , 1.0,
+ 1 'BLAMB1 ' , ' B color lambda 1 ' , 4300.0,
+ 2 'BLAMB2 ' , ' B color lambda 2 ' , 4700.0,
+ 3 'BLUM_A ' , ' frac. B lumin.for A' , 0.0,
+ 4 'BLUM_B ' , ' frac. B lumin.for B' , 0.0,
+ 5 'BLUM_C ' , ' frac. B lumin.for C' , 0.0,
+ 6 'BALB_A ' , ' albedo B for A ' , 1.0,
+ 7 'BALB_B ' , ' albedo B for B ' , 1.0,
+ 8 'BALB_C ' , ' albedo B for C ' , 1.0,
+ 1 'VLAMB1 ' , ' V color lambda 1 ' , 5300.0,
+ 2 'VLAMB2 ' , ' V color lambda 2 ' , 5700.0,
+ 3 'VLUM_A ' , ' frac. V lumin.for A' , 0.0,
+ 4 'VLUM_B ' , ' frac. V lumin.for B' , 0.0,
+ 5 'VLUM_C ' , ' frac. V lumin.for C' , 0.0,
+ 6 'VALB_A ' , ' albedo V for A ' , 1.0,
+ 7 'VALB_B ' , ' albedo V for B ' , 1.0,
+ 8 'VALB_C ' , ' albedo V for C ' , 1.0,
+ 1 'RLAMB1 ' , ' R color lambda 1 ' , 6500.0,
+ 2 'RLAMB2 ' , ' R color lambda 2 ' , 7500.0,
+ 3 'RLUM_A ' , ' frac. R lumin.for A' , 0.0,
+ 4 'RLUM_B ' , ' frac. R lumin.for B' , 0.0,
+ 5 'RLUM_C ' , ' frac. R lumin.for C' , 0.0,
+ 6 'RALB_A ' , ' albedo R for A ' , 1.0,
+ 7 'RALB_B ' , ' albedo R for B ' , 1.0,
+ 8 'RALB_C ' , ' albedo R for C ' , 1.0,
+ 1 'ILAMB1 ' , ' I color lambda 1 ' , 7500.0,
+ 2 'ILAMB2 ' , ' I color lambda 2 ' , 8500.0,
+ 3 'ILUM_A ' , ' frac. I lumin.for A' , 0.0,
+ 4 'ILUM_B ' , ' frac. I lumin.for B' , 0.0,
+ 5 'ILUM_C ' , ' frac. I lumin.for C' , 0.0,
+ 6 'IALB_A ' , ' albedo I for A ' , 1.0,
+ 7 'IALB_B ' , ' albedo I for B ' , 1.0,
+ 8 'IALB_C ' , ' albedo I for C ' , 1.0 /
+ DATA (NAMPAR(I),DESCPAR(I),PARDEF(I),I=101,116)
+ 1 /'ALBEDO ' , ' bolometric albedo ' , 1.0,
+ 2 'MAXITER ' , ' maximum iteration ' , 5.0,
+ 3 'PRECISION ' , ' convergency check ' , 0.001,
+ 4 'COARSE ' , ' coarsing factor ' , 1.0,
+ 1 'I ' , ' incl.of orbit plane' , 90.0,
+ 2 'MRATIO ' , ' mass ratio Mb/Ma ' , 1.0,
+ 3 'ECC ' , ' orbit eccentricity ' , 0.0,
+ 4 'PREC ' , ' Newton-Rapson prec.' , 1.E-7,
+ 5 'NULL ' , ' Unused ' , 0.0,
+ 6 'NULL ' , ' Unused ' , 0.0,
+ 1 'NUMPHASES ' , ' num.of equis.phases' , 8.0,
+ 2 'PHASEUNIT ' , ' unit for input phas' , 0.0,
+ 3 'NORM ' , ' Light curve norm. ' , 1.0,
+ 1 'NUMCELLS ' , ' num.of grid points ' , 100.0,
+ 2 'R4PREC ' , ' sin,cos precision ' , 1.E-5,
+ 3 'NULL ' , ' Unused ' , 0.0 /
+C
+ DATA (NAMPAR(I),DESCPAR(I),PARDEF(I),I=117,MAXPAR)
+ 1 /'AMOUNT ' , ' Amount of printing ' , 10.0,
+ 2 'SCREEN ' , ' Output on screen ' , 1.0,
+ 3 'UNIT ' , ' Output on this unit' , 11.0,
+ 4 'MAP ' , ' Map print phas.step' , 1.0,
+ 5 'REFL ' , ' Reflection details ' , 10.0,
+ 6 'LIGHT ' , ' Light curve details' , 4.0,
+ 1 'AMOUNT ' , ' Amount of printing ' , 10.0,
+ 2 'SCREEN ' , ' Output on screen ' , 1.0,
+ 3 'UNIT ' , ' Output on this unit' , 11.0,
+ 4 'MAP ' , ' Map plot phase step' , 1.0,
+ 5 'FILE ' , ' Output ASCII file ' , 12.0,
+ 6 'NULL ' , ' unused parameter ' , 0.0 /
+C
+C .......................................................
+ NPNTMX=MAXPT
+ NPNTMXC=MAXPTC
+ NALLMX=MAXALLIN
+ NFASIMX=MAXFASI
+ NBANDE=MAXBANDE
+ NPAR=MAXPAR
+ NFLAG=MAXFLAG
+ NTITLEMX=MAXTITLE
+ NFILT=MAXFILT
+ NLAM=MAXLAM
+C
+C .......................... elapsed time computation
+ TEMPO0=SECNDS(0.0) ! VAX VMS, NOT ANSI STANDARD
+ TEMPREC=0.0
+ TEMPIO=0.0
+ TEMPIOP=0.0
+ TEMPDRAW=0.0
+ TEMPREF=0.0
+ TEMPCOL=0.0
+ TEMPNORM=0.0
+ TEMPLUC=0.0
+C
+C .......................... Default setting
+ DO 3 I=1,MAXFLAG
+ 3 IFLAG(I)=IFLAGDEF(I)
+ DO 4 I=1,MAXPAR
+ 4 PAR(I)=PARDEF(I)
+C .......................... Colour band definitions
+ CALL FILFILT
+ DO 5 J=1,3
+ DO 5 I=1,MAXFILT
+ 5 COMPFRAC(I,J)=0.0
+C .......................... At the beginning there aren't stars
+ NTITLE=0
+ NPNT=0
+ NPNT1=0
+ NPNT2=0
+ NPNT3=0
+ N1I=0
+ N1F=0
+ N2I=0
+ N2F=0
+ N3I=0
+ N3F=0
+ NPNTC=0
+ NPNT1C=0
+ NPNT2C=0
+ NPNT3C=0
+ N1IC=0
+ N1FC=0
+ N2IC=0
+ N2FC=0
+ N3IC=0
+ N3FC=0
+C ...................... At the beginning there isn't project.grid
+ NCELL=0
+ DCELL=0
+C ...................... Testing some internal dimension
+ IF(NFILT.NE.NBANDE) THEN
+ WRITE(6,990) NFILT,NBANDE
+ 990 FORMAT(' WARNING! The number of colour bands in commons ',
+ 1 ' light and filtri is different:',2I5)
+ ENDIF
+C
+ WRITE(6,999)
+ 999 FORMAT(1H1,20X,' CLOSE BINARY SYSTEM ANALYSIS PROGRAM'/
+ 1 20X,' Version 0.0 ( debugging in progress )'//)
+C
+C ------------
+C =================================== => | INPUT |
+C ------------
+ 12 CONTINUE
+ CALL INPUT(NAMFLAG,DESCFLAG,IFLAGDEF,NAMPAR,DESCPAR,PARDEF,
+ 1 IPPAR1,IPPAR2,STOP)
+ IF(STOP) THEN
+C Elapsed time computation
+ IF(IFLAG(20).GT.0) THEN
+ TEMPO=SECNDS(TEMPO0)
+ IF(PRINT6)
+ 1 WRITE(6,8010) TEMPIO,TEMPIOP,TEMPDRAW,TEMPREF,
+ 2 TEMPCOL,TEMPNORM,TEMPLUC,TEMPO
+ 8010 FORMAT(///' Total elapsed time for I/O: ',G20.5/
+ 1 ' Total elapsed time for plotting: ',G20.5/
+ 2 ' Total elapsed time for surface drawing: ',G20.5/
+ 3 ' Total elapsed time for reflection : ',G20.5/
+ 4 ' Total elapsed time for colour band flux: ',G20.5/
+ 4 ' Total elapsed time for flux normalization:',G20.5/
+ 5 ' Total elapsed time for light to observer:',G20.5/
+ 6 ' Total elapsed time for this run: ',G20.5///)
+ IF(PRINTFILE)
+ 1 WRITE(NFILE,8010) TEMPIO,TEMPIOP,TEMPDRAW,TEMPREF,
+ 2 TEMPCOL,TEMPNORM,TEMPLUC,TEMPO
+ ENDIF
+ STOP
+ ENDIF
+C -----------------
+C =================================== => | RUN BEGINS |
+C -----------------
+ 300 CONTINUE ! 300: unreferenced label
+C MODFLAG is used to set off some flags after the corresponding comput.:
+C run flag: n.15 - input phases: n.12 -
+C stars and disk drawing: n.1,2,3
+C reflection and zero reflection source n.9,10
+ MODFLAG(15)=1
+ MODFLAG(12)=0
+ MODFLAG(1)=0
+ MODFLAG(2)=0
+ MODFLAG(3)=0
+ MODFLAG(9)=0
+ MODFLAG(10)=0
+C .......... sets print flags
+ IF(IFLAG(20).LE.0) THEN
+ NPRINT=0
+ NFILE=6
+ NPRINTR=0
+ ELSE
+ K1=IPPAR1(20)
+ AMOUNT=PAR(K1)
+C Screen output
+ NFILE6=INT(PAR(K1+1))
+ IF(NFILE6.GT.0) THEN
+ PRINT6=.TRUE.
+ ELSE
+ PRINT6=.FALSE.
+ ENDIF
+C Output on unit NFILE
+ NFILE=INT(PAR(K1+2))
+ IF(NFILE.GT.0.AND.NFILE.LE.99) THEN
+ PRINTFILE=.TRUE.
+ ELSE
+ PRINTFILE=.FALSE.
+ ENDIF
+C Phase step for printing maps
+ NPRINT=INT(PAR(K1+3))
+C Light curve details printing
+ NPRINTL=INT(PAR(K1+5))
+C Reflection hystory printing
+ NPRINTR=INT(PAR(K1+4))
+ IF(PRINTFILE) WRITE(NFILE,999)
+ ENDIF
+C .......... sets plot flags
+ IF(IFLAG(21).LE.0) THEN
+ NPLOT=0
+ NFILEP=6
+ NPRINTR=0
+ ELSE
+ K1=IPPAR1(21)
+ AMOUNTP=PAR(K1)
+C Screen output
+ NFILE6P=INT(PAR(K1+1))
+ IF(NFILE6P.GT.0) THEN
+ PLOT6=.TRUE.
+ ELSE
+ PLOT6=.FALSE.
+ ENDIF
+C Output on unit NFILEP
+ NFILEP=INT(PAR(K1+2))
+ IF(NFILEP.GT.0.AND.NFILEP.LE.99) THEN
+ PLOTFILE=.TRUE.
+ ELSE
+ PLOTFILE=.FALSE.
+ ENDIF
+C Phase step for printing maps
+ NPLOT=INT(PAR(K1+3))
+C Otput on ASCII file for external plotting
+ NASCIIP=INT(PAR(K1+4))
+ ENDIF
+C
+C ..................... sin.cos precision parameter
+C Sometimes sin,cos < r4prec is assumed 0, avoiding rounding error in
+C computing sin and cos near zero, without using double precision.
+ K1=IPPAR1(13)+1
+ R4PREC=PAR(K1)
+C ---------------------
+C .......................................... Input of phase values
+C ---------------------
+ IF(IFLAG(12).GT.0) THEN
+ K1=IPPAR1(12)
+ NFILEF=INT(PAR(K1+1))
+ IF(NFILEF.GT.0) THEN
+ CALL FASREAD(NFASI,NFILEF,NFILE,NFASIMX,FASE)
+ MODFLAG(12)=1
+ PAR(K1)=NFASI
+ ELSE
+ NFASI=INT(PAR(K1))
+ IF(NFASI.GT.0) THEN
+ IF(NFASI.GT.NFASIMX) THEN
+ WRITE(6,3000) NFASIMX,NFASIMX
+ IF(PRINTFILE)WRITE(NFILE,3000)NFASI,NFASIMX,NFASIMX
+ 3000 FORMAT(' ERROR! Too mutch phase values requested:',I5/
+ 1 ' It is set to the maximum value :',I5/
+ 2 ' If you need a greater value go into the FORTRAN list to',
+ 3 ' change: MAXFASI=',I5)
+ NFASI=NFASIMX
+ PAR(K1)=NFASI
+ ENDIF
+ CALL FASGRID(NFASI,FASE)
+ MODFLAG(12)=1
+ ENDIF
+ ENDIF
+ ENDIF
+C ...............................................................
+C
+C I/O Elapsed time computation
+ IF(IFLAG(20).GT.0) THEN
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPIO=TEMPIO+TEMPD
+ IF(PRINT6) WRITE(6,8001) TEMPD,TEMPIO
+C IF(PRINTFILE) WRITE(NFILE,8001) TEMPD,TEMPIO
+ 8001 FORMAT(/' Elapsed time for I/O:',12X,G20.5,5X,' Total:',G20.5/)
+ ENDIF
+C
+C ................ If Roche model computes Lagrange Points L1,L2,L3
+ IF(PARS(1,1).EQ.2..OR.PARS(2,1).EQ.2..OR.PARS(3,1).EQ.2.) THEN
+ K1=IPPAR1(11)
+ Q=PAR(K1+1)
+C IF(Q.GT.1.) THEN ! The following isn't true
+C Q1=1./Q
+C WRITE(6,3020)Q,Q1
+C IF(PRINTFILE) WRITE(NFILE,3020)Q,Q1
+C 3020 FORMAT(/' ERROR ! For Roche model object A must be the'
+C 1 ' most massive'/' q:',G14.7,' changed to:',G14.7)
+C Q=Q1
+C PAR(K1+1)=Q
+C ENDIF
+ PRECISION=PAR(IPPAR1(11)+3)
+ IF(PRECISION.LT.1.E-7) WRITE(6,3021) PRECISION
+ 3021 FORMAT(//' WARNING ! Requested Newton-Rapson precision ',
+ 1 ' may be too small to be used with real*4 numbers:',G14.7//)
+ CALL LAGR(Q,PRECISION)
+ IF(PRINT6) WRITE(6,3022) AL1,OMEGAL1,AL2,OMEGAL2,AL3,OMEGAL3
+ 3022 FORMAT(/' Computed Lagrange points and potential for',
+ 1 ' critical lobes:'/,5X,'L1=',G14.7,' omega=',G14.7/
+ 2 5X,'L2=',G14.7,' omega=',G14.7/
+ 3 5X,'L3=',G14.7,' omega=',G14.7/)
+ ENDIF
+C -------------------
+C ========== => | SURFACE DRAWING |
+C -------------------
+C
+C ...............................................................
+C | The surface of each object (the 2 stars and the disk), is |
+C | represented by a number of "surface elements", which are |
+C | stored in common /surf/ as a list containing, for each |
+C | element, the area, the temperature, the flux etc. |
+C ...............................................................
+C
+C ------------------------------------------loop on stars and disk
+ DO 30 IS=1,3
+C If the object number IS is not to be drawn goes on
+ IF(IFLAG(IS).LE.0) GOTO 30
+ MODFLAG(IS)=1
+C The reflection has to be recomputed with new grid
+ IF(IFLAG(9).LT.0) IFLAG(9)=1
+ IF(IFLAG(10).LT.0) IFLAG(10)=1
+C
+C If the object had been drawn before,
+C deletes the old space, rearranging data.
+C ALL data in common/surf/are mouved,also if not jet computed.
+C Also /surfc/ is compressed.
+ IF(NP123(IS).GT.0.AND.NIF(1,IS).GT.0.AND.NIF(2,IS).GT.0)
+ 1 CALL CANCEL(IS,NFILE)
+C
+ NPS=NP123(IS)
+ NIS=NIF(1,IS)
+ NFS=NIF(2,IS)
+ IF(NIS.LE.0) NIS=NPNT+1
+ NPSC=NP123C(IS)
+ NISC=NIFC(1,IS)
+ NFSC=NIFC(2,IS)
+ IF(NISC.LE.0) NISC=NPNTC+1
+C
+C ......... for sphere or Roche model test and set number of intervals
+C
+ IF(PARS(1,IS).EQ.1..OR.PARS(1,IS).EQ.2.) THEN
+C mesh parameter numtheta: number of fine mesh division.
+ NUMTHETA=INT(PARS(7,IS))
+C Tests if numtheta is consistent, it must be odd > 3
+ IF(NUMTHETA.LT.4) NUMTHETA=5
+ IF(MOD(NUMTHETA,2).EQ.0) NUMTHETA=NUMTHETA+1
+C Tests if the coarsing factor is consistent, it must be odd
+C and a submultiple of numtheta-1 , if not so nthf and
+C and numtheta are changed.
+ K1=IPPAR1(9)+3
+ 301 NTHF=INT(PAR(K1))
+C Makes NTHF odd > 0
+ IF(NTHF.EQ.0) NTHF=1
+ IF(MOD(NTHF,2).EQ.0) NTHF=NTHF+1
+C Makes numtheta consistent with NTHF
+ 31 CONTINUE
+ NTHRESTO=MOD(NUMTHETA-1,NTHF)
+C This is equal to: MOD(NUMTHETA+2*INT(NTHF/2),NTHF)
+ IF(NTHRESTO.NE.0) THEN
+C changes numtheta
+ NUMTHETA=NUMTHETA + NTHF - NTHRESTO
+C If not so, numtheta is made odd (by adding nthf, odd):
+ IF(MOD(NUMTHETA,2).EQ.0) NUMTHETA=NUMTHETA+NTHF
+C If numtheta too great
+ IF(NUMTHETA.GT.NPNTMX) THEN
+ IF(NUMTHETA-NTHF.GT.5) THEN
+C decrease numtheta
+ NUMTHETA=NUMTHETA-NTHF
+ ELSE
+C NTHF is decreased by 2 (it's odd)
+ NTHF=NTHF-2
+ GOTO 31
+ ENDIF
+ ENDIF
+ ENDIF
+C Warnng message and changes input coarsing factor
+ IF(INT(PAR(K1)).NE.NTHF) THEN
+ WRITE(6,3053) PAR(K1),NTHF,IS
+ IF(PRINTFILE) WRITE(NFILE,3053) PAR(K1),NTHF,IS
+ 3053 FORMAT(/' WARNING! Coarsing factor:',G15.7,
+ 1 ' has been changed to:',I6/' The coarsing factor',
+ 2 ' must be odd > 0 and consistent with the mesh',
+ 3 ' parameter.',5X,'Object:',I3)
+C PAR(K1)=NTHF lascio il valore precedente, mentre usa nthf in sfera
+ ENDIF
+C Warning message and changes input numtheta
+ IF(INT(PARS(7,IS)).NE.NUMTHETA) THEN
+ WRITE(6,3054) PARS(7,IS),NUMTHETA,NTHF,IS
+ IF(PRINTFILE) WRITE(NFILE,3054)
+ 1 PARS(7,IS),NUMTHETA,NTHF,IS
+ 3054 FORMAT(/' WARNING! MESH parameter:',G15.7,
+ 1 ' changed to:',I6/
+ 1 ' For spheres and Roche it must be odd > 4',
+ 1 ' and MESH-1 multiple of coarsing factor:',I5,
+ 2 3X,'Object:',I3)
+ PARS(7,IS)=NUMTHETA
+ ENDIF
+ ELSE IF(PARS(1,IS).EQ.0.0) THEN
+C ........................... similar tests for rectangle
+ NUMTHETA=INT(PARS(7,IS))
+ IF(NUMTHETA.LT.0) NUMTHETA=1
+ K1=IPPAR1(9)+3
+ NTHF=INT(PAR(K1))
+C Makes NTHF odd > 0
+ IF(NTHF.EQ.0) NTHF=1
+C The coarsing factor must be a submultiple of theta
+ NTHRESTO=MOD(NUMTHETA,NTHF)
+ IF(NTHRESTO.NE.0) NUMTHETA=NUMTHETA+NTHF-NTHRESTO
+ IF(NUMTHETA.NE.INT(PARS(7,IS)) ) THEN
+ WRITE(6,3062) PARS(7,IS),NUMTHETA,IS
+ IF(PRINTFILE) WRITE(NFILE,3062) PARS(7,IS),NUMTHETA,IS
+ PARS(7,IS)=NUMTHETA
+ ENDIF
+ ELSE IF(PARS(1,IS).EQ.3.0) THEN
+C ..................... number of intervals for disk
+ NUMTHETA=INT(PARS(7,IS))
+ IF(NUMTHETA.LT.0) NUMTHETA=1
+ K1=IPPAR1(9)+3
+ NTHF=INT(PAR(K1))
+C Makes NTHF odd > 0
+ IF(NTHF.EQ.0) NTHF=1
+C IF(MOD(NTHF,2).EQ.0) THEN
+C PAR(K1)=PAR(K1)+1
+C WRITE(6,3052) NTHF,PAR(K1),IS
+C IF(PRINTFILE) WRITE(NFILE,3052) NTHF,PAR(K1),IS
+C NTHF=NTHF+1
+C3052 FORMAT(/' WARNING! The input coarsing factor:',
+C 1 I5,' has been changed to the odd value:',G15.7,
+C 2 2X,'Object:',I3)
+C ENDIF
+C The coarsing factor must be a submultiple of theta
+ NTHRESTO=MOD(NUMTHETA,NTHF)
+ IF(NTHRESTO.NE.0) NUMTHETA=NUMTHETA+NTHF-NTHRESTO
+ IF(NUMTHETA.NE.INT(PARS(7,IS)) ) THEN
+ WRITE(6,3062) PARS(7,IS),NUMTHETA,IS
+ IF(PRINTFILE) WRITE(NFILE,3062) PARS(7,IS),NUMTHETA,IS
+ 3062 FORMAT(/' WARNING! The mesh parameter :',
+ 1 G15.7,' has been changed to the value:',I5/
+ 2 ' Consistent with the coarsing factor',5X,'Object:',I3)
+ PARS(7,IS)=NUMTHETA
+ ENDIF
+ ENDIF
+C ------------
+C ............. SINGLE POINT ,in a rectangular interval
+C ------------
+ IF(PARS(1,IS).EQ.0.0) THEN
+C
+C
+C Test other values: radius: i.e. X half-edge
+ RAGGIO=PARS(2,IS)
+ IF(RAGGIO.LT.0.0) THEN
+ RAGGIO=0.
+ WRITE(6,3010) IS,PARS(2,IS),RAGGIO
+ IF(PRINTFILE) WRITE(NFILE,3010) IS,PARS(2,IS),RAGGIO
+ 3010 FORMAT(' ERROR! : Object:',I3,' Absurd radius value:',G15.7,
+ 1 ' changed to:',G15.7)
+ PARS(2,IS)=0.0
+ ENDIF
+C
+C Y half-edge
+ RAGGIOI=PARS(15,IS)
+ IF(RAGGIOI.LT.RAGGIO) THEN
+ RAGGIOI=RAGGIO
+ WRITE(6,3012) IS,PARS(15,IS),RAGGIOI
+ IF(PRINTFILE) WRITE(NFILE,3012) IS,PARS(15,IS),RAGGIOI
+ 3012 FORMAT(' Object:',I3,' Y Input radius ( < X radius):',G15.7,
+ 1 ' changed to X radius:',G15.7)
+ ELSE IF(MOD(RAGGIOI,RAGGIO).GT.R4PREC) THEN
+ RAGGIOI=RAGGIO*INT(RAGGIOI/RAGGIO)
+ WRITE(6,3014) IS,PARS(15,IS),RAGGIOI
+ IF(PRINTFILE) WRITE(NFILE,3014) IS,PARS(15,IS),RAGGIOI
+ 3014 FORMAT(' WARNING! Object:',I3/' Y Input radius:',G15.7,
+ 1 ' changed to :',G15.7,' multiple of X radius')
+ ENDIF
+C
+C NPSSC is the number of the first void coarse, it is updated
+C by CYL, NPSS is computed by CYL as the number of fine made by CYL
+ NPSSC=NISC
+ NPNTMX1=NPNTMX-NIS+1
+ ALZO =(PI/180.)*PARS(18,IS)
+ ROTAZ=(PI/180.)*PARS(19,IS)
+ PKRI=PARS(20,IS)
+C
+ CALL POINT(PRINT6,PRINTFILE,NFILE,NUMTHETA,NTHF,PKRI,IS,
+ 1 NPSS,NPSSC,RAGGIO,RAGGIOI,
+ 2 PARS(3,IS),PARS(4,IS),PARS(5,IS),
+ 3 ALZO,ROTAZ,R4PREC,NPNTMX1,X(NIS),Y(NIS),Z(NIS),
+ 4 G(NIS),GX(NIS),GY(NIS),GZ(NIS),
+ 5 FX(NIS),FY(NIS),FZ(NIS),TX(NIS),TY(NIS),TZ(NIS),
+ 6 AREA(IS),A(NIS),NUMCOARSE(NIS),NUMOBJ(NIS))
+C
+C NUMOBJ is assigned into subroutine POINT; not outside, as in SFERA2.
+C
+C -----------------
+C ............. SPHERICAL SURFACE, drawn by equi-area intervals
+C -----------------
+ ELSE IF( PARS(1,IS).EQ.1.0) THEN
+C
+C Test other values
+ RAGGIO=PARS(2,IS)
+ IF(RAGGIO.LT.0) THEN
+ RAGGIO=0.
+ WRITE(6,3010) IS,PARS(2,IS),RAGGIO
+ IF(PRINTFILE) WRITE(NFILE,3010) IS,PARS(2,IS),RAGGIO
+ PARS(2,IS)=0.0
+ ENDIF
+C
+C NPSSC is updated and NPSS is computed by sfera2:
+C if you are drawing surfaces for the first time, then NPSSC
+C is the next free coarse surf element (total number of coarse points
+C (for all objects) +1 ).Otherwise the first coarse of next object.
+C NPSS is the number of fine surf.interval for this object.
+C NCORD is the cord option, NPNTMX1 the available space for surf.el.
+ NPSSC=NISC
+ NCORD=INT(PARS(14,IS))
+ NPNTMX1=NPNTMX-NIS+1
+C
+ CALL SFERA2(PRINT6,PRINTFILE,NFILE,NUMTHETA,NTHF,NCORD,NPSS,
+ 1 NPSSC,RAGGIO,PARS(3,IS),PARS(4,IS),PARS(5,IS),
+ 2 NPNTMX1,X(NIS),Y(NIS),Z(NIS),G(NIS),GX(NIS),
+ 3 GY(NIS),GZ(NIS),FX(NIS),FY(NIS),FZ(NIS),
+ 4 TX(NIS),TY(NIS),TZ(NIS),
+ 5 AREA(IS),A(NIS),NUMCOARSE(NIS),AUS,AUS1 )
+C
+C NUMOBJ, number of object for each surf. el., is used in LUCE1
+C to distinguish between the surf.el. of different objects and printings
+ CALL FILL(NPSS,IS,NUMOBJ(NIS))
+C
+ ELSE IF (PARS(1,IS).EQ.2.) THEN
+C ----------
+C ................................... ROCHE LOBE drawing
+C ----------
+C The star A (main body) is drawn in (0,0,0) and
+C the star B in (1,0,0) (The distance between the stars
+C acts as the distance unit measure). The stars can be
+C shifted to an input given x point: PARS(3,is), this
+C is equivalent to change Omega, and will be faster.
+C This option is't valid for object 3
+ IF(IS.EQ.3) THEN
+ WRITE(6,3073)
+ 3073 FORMAT(' ERROR! Roche model not allowed for object 3,',
+ 1 ' object not drawn!')
+ GOTO 30
+ ENDIF
+C .....Testing y and z for the body center
+ IF(PARS(4,IS).NE.0.0.OR.PARS(5,IS).NE.0.) THEN
+ WRITE(6,3074) IS
+ 3074 FORMAT(' WARNING! ,for Roche model the object:',I3,
+ 1 ' must be in (X,0,0)')
+ PARS(4,IS)=0.0
+ PARS(5,IS)=0.0
+ ENDIF
+C
+ Q=PAR(IPPAR1(11)+1)
+ PRECISION=PAR(IPPAR1(11)+3)
+ IF(PRECISION.LT.1.E-7) WRITE(6,3021) PRECISION
+C ... test lobes-x axis intersection, computes omega from R if not given
+ CALL LOBES(Q,PRECISION,PARS(6,IS),PARS(2,IS),IS)
+ XSHIFT=PARS(3,IS)
+C for object 2 the star is drawn n (1,0,0) : sets the correct shift:
+ IF(IS.EQ.2) XSHIFT=XSHIFT-1.
+ OMEGA=PARS(6,IS)
+ NPSSC=NISC
+ NCORD=INT(PARS(14,IS))
+ NPNTMX1=NPNTMX-NIS+1
+ IF(IS.EQ.1) THEN
+C the first point XA is near Lagr.point L1
+ XA=XA1
+ XB=XA2
+ RIS=0.0
+ IF(PRINT6) WRITE(6,3080) IS,XA1,XA2,XA
+ 3080 FORMAT(/' Surface-X axis intersection for object:',I2,
+ 1 ' computed as:',G14.7,2X,G14.7/
+ 2 ' Surface drawing starts from point:',G14.7/)
+C
+ CALL ROCHES(PRINT6,PRINTFILE,NFILE,NUMTHETA,NTHF,NCORD,NPSS,
+ 1 NPSSC,XA,XB,XSHIFT,OMEGA,Q,RIS,PRECISION,
+ 2 NPNTMX1,X(NIS),Y(NIS),Z(NIS),G(NIS),GX(NIS),
+ 3 GY(NIS),GZ(NIS),FX(NIS),FY(NIS),FZ(NIS),
+ 4 TX(NIS),TY(NIS),TZ(NIS),
+ 5 AREA(IS),A(NIS),NUMCOARSE(NIS),AUS,AUS1,ROCHESF1)
+ ELSE
+C the last point XB is near Lagr.point L1
+ XA=XB2
+ XB=XB1
+ RIS=1.0
+ IF(PRINT6) WRITE(6,3080) IS,XB1,XB2,XA
+C
+ CALL ROCHES(PRINT6,PRINTFILE,NFILE,NUMTHETA,NTHF,NCORD,NPSS,
+ 1 NPSSC,XA,XB,XSHIFT,OMEGA,Q,RIS,PRECISION,
+ 2 NPNTMX1,X(NIS),Y(NIS),Z(NIS),G(NIS),GX(NIS),
+ 3 GY(NIS),GZ(NIS),FX(NIS),FY(NIS),FZ(NIS),
+ 4 TX(NIS),TY(NIS),TZ(NIS),
+ 5 AREA(IS),A(NIS),NUMCOARSE(NIS),AUS,AUS1,ROCHESF2)
+ ENDIF
+C Assign the object number to the surface elements
+ CALL FILL(NPSS,IS,NUMOBJ(NIS))
+C
+C ----
+C .................................. DISK drawing
+C ----
+ ELSE IF (PARS(1,IS).EQ.3.) THEN
+C
+C Test other values
+ RAGGIO=PARS(2,IS)
+ IF(RAGGIO.LT.0) THEN
+ RAGGIO=0.
+ WRITE(6,3010) IS,PARS(2,IS),RAGGIO
+ IF(PRINTFILE) WRITE(NFILE,3010) IS,PARS(2,IS),RAGGIO
+ PARS(2,IS)=0.0
+ ENDIF
+ RAGGIOI=PARS(15,IS)
+ IF(RAGGIOI.LT.0.0) THEN
+ RAGGIOI=0.0
+ WRITE(6,3010) IS,PARS(15,IS),RAGGIOI
+ IF(PRINTFILE) WRITE(NFILE,3010) IS,PARS(15,IS),RAGGIOI
+ PARS(2,IS)=0.0
+ ENDIF
+C
+C NPSSC is the number of the first void coarse, it is updated
+C by CYL, NPSS is computed by CYL as the number of fine made by CYL
+ NPSSC=NISC
+ NCORD=INT(PARS(14,IS))
+ NPNTMX1=NPNTMX-NIS+1
+ ALZO =(PI/180.)*PARS(18,IS)
+ ROTAZ=(PI/180.)*PARS(19,IS)
+ PKRI=PARS(20,IS)
+C
+ CALL CYL(PRINT6,PRINTFILE,NFILE,NUMTHETA,NTHF,NCORD,PKRI,IS,
+ 1 NPSS,NPSSC,RAGGIO,RAGGIOI,PARS(16,IS),
+ 2 PARS(17,IS),PARS(3,IS),PARS(4,IS),PARS(5,IS),
+ 3 ALZO,ROTAZ,R4PREC,NPNTMX1,X(NIS),Y(NIS),Z(NIS),
+ 4 G(NIS),GX(NIS),GY(NIS),GZ(NIS),
+ 5 FX(NIS),FY(NIS),FZ(NIS),TX(NIS),TY(NIS),TZ(NIS),
+ 6 AREA(IS),A(NIS),NUMCOARSE(NIS),NUMOBJ(NIS),
+ 7 AUS,AUS(NTHF+1),AUS(2*NTHF+1),AUS(3*NTHF+1),
+ 8 AUS(4*NTHF+1) )
+C
+C NUMOBJ is assigned into subroutine CYL; not outside, as in SFERA2.
+C
+ ENDIF
+C Assign the parameters of the surface element list
+C for the drawn object
+ NIF(1,IS)=NIS
+ NIF(2,IS)=NIS+NPSS-1
+ NP123(IS)=NPSS
+ NPNT=NPNT+NPSS
+C nifc(2,is)=0 means coarse map done, but NOT coarse surf elem.
+ NIFC(1,IS)=NISC
+ NIFC(2,IS)=0
+ NP123C(IS)=NPSSC-NISC
+ NPNTC=NPNTC+NP123C(IS)
+C
+C -----------------------------------------------------
+C ................GRAVITY DARKENED TEMPERATURE PROFILE AND SURF.EL FLUX
+C -----------------------------------------------------
+C ( Sets FB(i),T(i),i=nis,nfs)
+C
+C If g , gravity darkening parameter, is set
+ IF(PARS(9,IS).NE.0.0) THEN
+C Computes gravity darkened temperature given by beta:
+C Pars(9,is) is the g gravity dark. coeff. definition by Wilson
+C different from Lucey and Von Zeipel beta coefficient
+C pars(10,is) is the input pole temperature= max T
+ BETA=PARS(9,IS)*0.25
+ CALL GRAVITY1(KK,NPSS,BETA,PARS(10,IS),G(NIS),T(NIS))
+C Store pole values for printing
+ KK=NIS+KK-1
+ XPOLE(IS)=X(KK)
+ YPOLE(IS)=Y(KK)
+ ZPOLE(IS)=Z(KK)
+ TPOLE(IS)=T(KK)
+ RPOLE(IS)=SQRT((X(KK)-PARS(3,IS))**2 +(Y(KK)-PARS(4,IS))**2+
+ 1 (Z(KK)-PARS(5,IS))**2 )
+ GXPOLE(IS)=GX(KK)
+ GYPOLE(IS)=GY(KK)
+ GZPOLE(IS)=GZ(KK)
+C
+C and bolometric flux is given by T**4 in each point
+C flux = caT**4/(4PI)=erg/cm**2/sec/sterad
+C to shorten calculation we could assume here c*a*T**4/(4*pi)=1
+ CALL LBOLT(NPSS,FB(NIS),T(NIS),A(NIS))
+ ELSE
+C uniform temperature and flux in each surface element
+ EFFE=PARS(10,IS)**4 *AC4PI
+ CALL FILL(NPSS,PARS(10,IS),T(NIS))
+ CALL FILL1(NPSS,EFFE,FB(NIS),A(NIS))
+ ENDIF
+C --------------------------------
+C ....................... LIMB DARKENING AND ALBEDO COEFF.
+C --------------------------------
+C (Sets XLIMB(i),ALB(i),i=nis,nfs)
+C eventually a T-dependence could be inserted here
+ DARKLIM=PARS(13,IS)
+ ALBIS=PARS(12,IS)
+ CALL LIMB(NPSS,DARKLIM,ALBIS,T(NIS),XLIMB(NIS),ALB(NIS))
+C
+C Drawing elapsed time computation
+ IF(IFLAG(20).GT.0) THEN
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPDRAW=TEMPDRAW+TEMPD
+ IF(PRINT6) WRITE(6,8002) TEMPD,TEMPDRAW
+C IF(PRINTFILE) WRITE(NFILE,8002) TEMPD,TEMPDRAW
+ 8002 FORMAT(/' Elapsed time for surface drawing:',
+ 1 G20.5,5X,' Total:',G20.5/)
+ ENDIF
+C
+ 30 CONTINUE
+C ---------------------------------------- end of loop on stars
+C
+C
+C
+C --------------------------
+C ============== => | REFLECTION |
+C --------------------------
+ IF(IFLAG(9).GT.0) THEN
+ MODFLAG(9)=1
+C In this case only obtain frifl from friflc from previous run
+ IF(IFLAG(9).EQ.2) GOTO 320
+C
+C Look if, for some object, the coarse map must be used
+C to compute the coarse surf.elements in common /surfc/:
+ DO 32 IS=1,3
+ IF(NIFC(2,IS).LE.0.AND.NP123C(IS).GT.0) THEN
+ IF(NPNTC.GT.NPNTMXC) THEN
+ WRITE(6,3100) IS,NP123C(IS),NPNTMXC
+ IF( PRINTFILE) WRITE(NFILE,3100) IS,NP123C(IS),NPNTMXC
+ 3100 FORMAT(' ERROR: for object',I2,' I can''t use',I6,
+ 1 ' coarse surface elements for reflection.'/
+ 2 ' Reduce the number of coarse surface elements for some object,'/
+ 3 ' or go into the FORTRAN list to change parameter MAXPTC, now:',
+ 4 I6)
+ ENDIF
+ CALL COARSE(PRINT6,PRINTFILE,NFILE,IS)
+ ENDIF
+ 32 CONTINUE
+C
+C Bolometric reflection is used to compute a bolometric
+C reflected flux for each surface element : FRIFL(NPNT)
+C An albedo parameter (ALBEDO) which wheights the reflection
+C efficiency is used for all the objects; a similar input
+C parameter exists for each object.
+ ALBEDO=PAR(IPPAR1(9))
+ MAXITER=INT(PAR(IPPAR1(9)+1))
+ PRECISION=PAR(IPPAR1(9)+2)
+C Computing reflected flux for each coarse surface elememt
+ CALL REFLECT(MAXITER,PRECISION,R4PREC,ALBEDO,NPRINTR,NFILE,
+ 1 AUS)
+C Computing total bolometric reflected lum. for each object
+ 320 CALL SUMREFL(NPNTC,NIFC,FBOLREF,FRIFLC)
+C From coarse to fine surf.el.reflected flux(FRIFLC=>FRIFL)
+ CALL FINEFRIFL
+C if black body is used for any color or object to obtain
+C flux fraction from t for each surf. element
+C T must changed in each surf. element after reflection
+ DO 33 IS=1,3
+ IF(IFLAG(IS).LE.0) GOTO 33
+ DO 33 IB=1,NBANDE
+ IF(IFLAG(IB+3).LE.0) GOTO 33
+ K1=IPPAR1(IB+3)+1
+ ALFRAC=PAR(K1+IS)
+ IF(ALFRAC.GE.0.0) GOTO 33
+C in this case T of single surf element is used (see loop 44 below)
+ CALL TCHANGE(NPNT,FB,FRIFL,T,A)
+ GOTO 330
+ 33 CONTINUE
+ 330 CONTINUE
+C If a T dependent limb darkening is used you should
+C iterate correcting t and then returninig to compute
+C reflection with the new limb.dark. coefficient.
+C This is not implemented being a too heavy computation.
+C
+C Reflection elapsed time computation
+ IF(PRINT6) THEN
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPREF=TEMPREF+TEMPD
+ IF (PRINT6) WRITE(6,8003) TEMPD,TEMPREF
+C IF (PRINTFILE) WRITE(NFILE,8003) TEMPD,TEMPREF
+ 8003 FORMAT(/' Elapsed time for reflection:',5X,G20.5,
+ 1 5X,' Total:',G20.5/)
+ ENDIF
+C
+ ENDIF
+C ........... Set reflection source to zero(Fbolref,frifl=0)
+ IF(IFLAG(10).GT.0) THEN
+ MODFLAG(10)=1
+ CALL NOREFL(NPNT,NIF,FBOLREF,FRIFL)
+ ENDIF
+C
+C ---------------------
+C .......................... TOTAL BOLOMETRIC FLUX
+C ---------------------
+C
+C Computing FBOL : total bolometric flux for
+C each object.
+C FBOLTOT, flux normalization factor is set =1.
+C Flux normalization is performad after colour band computation
+C
+ DO 35 IS=1,3
+ NIS=NIF(1,IS)
+ IF(NIS.LE.0) GOTO 35
+ NPSS=NP123(IS)
+ IF(NPSS.LE.0) GOTO 35
+ CALL TOTALE(NPSS,FBOL(IS),FB(NIS))
+ FBOLTOT(IS)=1.
+ 35 CONTINUE
+C
+C --------------------------
+C ================ => | SETS COLOUR BAND FLUX |
+C --------------------------
+C
+C ------------------------------------------loop on stars and disk
+ DO 40 IS=1,3
+ NIS=NIF(1,IS)
+ IF(NIS.LE.0) GOTO 40
+ NPSS=NP123(IS)
+ IF(NPSS.LE.0) GOTO 40
+C ............................ Sets the flux of each color band:
+C F(surf.el,colour)=FB(surf.el) * flux fraction in the colour
+C We have the following choices:
+C a)- alfrac>0 alfrac is used, flux is not computed
+C indepentent on alam1=alam2, alam1<>alam2
+C b)- alfrac<0 band flux is computed for each surf. element
+C T of single surf. el. is adjusted after reflection
+C in loop 33
+C c)- alfrac=0 colour band flux is computed from avarage T
+C (not corrected for reflection)
+C
+C 1 - monocromatic flux :alam1=alam2
+C 2 - band flux by integrating planK*colour filter
+C -----------------------------------------loop on colors
+ DO 44 IB=1,NBANDE
+C If the color band IB isn't requested goes on
+ IF(IFLAG(IB+3).LE.0) GOTO 44
+ K1=IPPAR1(IB+3)+1
+ ALFRAC=PAR(K1+IS)
+ TEMPER=PARS(10,IS)
+ ALAM1=PAR(IPPAR1(IB+3))
+ ALAM2=PAR(IPPAR1(IB+3)+1)
+ IF(ALFRAC.GT.0.0) THEN
+C input given flux fraction in the band F=(FB+FRIFL)*ALFRAC
+ CALL BANDFIL(NPSS,ALFRAC,FB(NIS),FRIFL(NIS),F(NIS,IB))
+ ELSE IF(ALAM1.EQ.ALAM2) THEN
+C monocromatic light curve
+ IF(ALFRAC.EQ.0.0) THEN
+C monocromatic flux FRACTION from pole temperature
+C plankl=plank*deltal with deltal=1.E-8cm
+ BANDL=PLANKL(ALAM1,TEMPER)/(AC4PI*TEMPER**4)
+C fills F with (FB+FRIFL)*bandl (FB=flux*area)
+ CALL BANDFIL(NPSS,BANDL,FB(NIS),FRIFL(NIS),F(NIS,IB))
+ ELSE IF(ALFRAC.LT.0.0) THEN
+C monocromatic flux fraction from t of each surf el.
+ CALL BANDATM(NPSS,ALAM1,T(NIS),F(NIS,IB),FB(NIS),
+ 1 FRIFL(NIS))
+C bandl=sum(f*area) tot flux in band for object
+ ENDIF
+ ELSE
+C if (alam1.ne.alam2) : integrating flux over colour bands
+ IF(ALFRAC.EQ.0.0) THEN
+C flux frac.by integrating plank( T pole) over the colour bands
+ CALL BANDAC(TOTBAND,TEMPER,NLAMBDA(IB),DELTALAM(IB),
+ 1 ALAMBDA(1,IB),TRASMISS(1,IB) )
+ BANDL=TOTBAND/(AC4PI*TEMPER**4)
+ CALL BANDFIL(NPSS,BANDL,FB(NIS),FRIFL(NIS),F(NIS,IB))
+ ELSE IF(ALFRAC.LT.0.0) THEN
+C flux frac.by integrating plank( T for each surf.el.)over bands
+C this can be a very heavy computation, don't use.
+ CALL BANDAT(NPSS,T(NIS),F(NIS,IB),FB(NIS),FRIFL(NIS),
+ 1 NLAMBDA(IB),DELTALAM(IB),ALAMBDA(1,IB),TRASMISS(1,IB))
+C bandl=sum(f*A) flux in band for object
+ ENDIF
+ ENDIF
+ CALL TOTALE(NPSS,TOTBAND,F(NIS,IB))
+ COMPFRAC(IB,IS)=TOTBAND
+C
+C Colour band flux elapsed time computation
+ IF(IFLAG(20).GT.0) THEN
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPCOL=TEMPCOL+TEMPD
+ IF(PRINT6) WRITE(6,8004) TEMPD,TEMPCOL
+C IF(PRINTFILE) WRITE(NFILE,8004)TEMPD,TEMPCOL
+ 8004 FORMAT(/' Elapsed time for color band:',5X,G20.5,5X,
+ 1 ' Total:',G20.5/)
+ ENDIF
+C
+ 44 CONTINUE
+C -------------------------------------- END of color loop
+ 40 CONTINUE
+C -------------------------------------- end of stars loop
+C
+C ---------------------
+C .......................... FLUX NORMALIZATION
+C ---------------------
+C
+C If a flux is input given for the object, normalizes the
+C BOL.flux + REFLECTION to the input flux.
+C or to BOL.flux (without reflection)
+C Reflection source is normalize in the same way. The colour
+C band flux is also normalized in a consistent way.
+C
+ DO 45 IS=1,3
+ NIS=NIF(1,IS)
+ IF(NIS.LE.0) GOTO 45
+ NPSS=NP123(IS)
+ IF(NPSS.LE.0) GOTO 45
+C .................... if normalization requested by input
+ IF(PARS(11,IS).LE.0.0) GOTO 45
+C
+C FACTOR=PARS(11,IS)/(FBOL(IS)+FBOLREF(IS))
+ FACTOR=PARS(11,IS)/FBOL(IS)
+C ..............bolometric flux normalization
+ CALL NORM(NPSS,FACTOR,FB(NIS))
+ FBOL(IS)=PARS(11,IS)
+ FBOLTOT(IS)=1./FACTOR
+C ..............reflected flux normalization ( if reflection is used )
+ IF(IFLAG(9).GT.0.AND.IFLAG(10).LE.0)
+ 1 CALL NORM(NPSS,FACTOR,FRIFL(NIS))
+ FBOLREF(IS)=FBOLREF(IS)*FACTOR
+C ..............color flux normalization ( only for used color bands)
+ DO 47 IB=1,NBANDE
+ IF(IFLAG(IB+3).LE.0) GOTO 47
+ CALL NORM(NPSS,FACTOR,F(NIS,IB))
+ COMPFRAC(IB,IS)=COMPFRAC(IB,IS)*FACTOR
+ 47 CONTINUE
+ 45 CONTINUE
+C
+C Normalization elapsed time computation
+ IF(IFLAG(20).GT.0) THEN
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPNORM=TEMPNORM+TEMPD
+ IF (PRINT6) WRITE(6,8015) TEMPD,TEMPNORM
+C IF (PRINTFILE) WRITE(NFILE,8015) TEMPD,TEMPNORM
+ 8015 FORMAT(/' Elapsed time for flux normalization:',5X,G20.5,
+ 1 5X,' Total:',G20.5/)
+ ENDIF
+C -----------------------------------
+C ================== => | COMPUTING LIGHT TO THE OBSERVER |
+C -----------------------------------
+C
+C ..if an object has been redrawn looks for new max and min x-y-z values
+ KTEST=.FALSE.
+ IF(IFLAG(1).GT.0.OR.IFLAG(2).GT.0.OR.IFLAG(3).GT.0) THEN
+C Maximum absolute value for x,y and z, which is the maximum
+C radius for the drawn system.
+ CALL MAXABS(CMAX1,NPNT,X)
+ CALL MAXABS(CMAX2,NPNT,Y)
+ CALL MAXABS(CMAX3,NPNT,Z)
+C The maximum possible linear dimension is the diagonal of x,y,z box
+ CELLMAX=SQRT(CMAX1*CMAX1+CMAX2*CMAX2+CMAX3*CMAX3)
+ CELLMIN=-CELLMAX
+ KTEST=.TRUE.
+ ENDIF
+C
+C ........................... Set grid points
+ K1=IPPAR1(13)
+ NGP=PAR(K1)
+ IF (NGP.LE.0.OR.NGP.GT.MAXCELL) THEN
+ NGP=MAXCELL
+ WRITE(6,3200) PAR(K1),NGP,NGP
+ IF(PRINTFILE) WRITE(NFILE,3200) PAR(K1),NGP,NGP
+ 3200 FORMAT(/' WARNING! The input number of grid points:',G15.7,
+ 1 ' can''t be greater than:',I5/
+ 2 ' If you need a greater value go to the FORTRAN list to',
+ 3 ' change MAXCELL=',I5)
+ ENDIF
+ IF (NCELL.NE.NGP.OR.KTEST) THEN
+ NCELL=NGP
+ DCELL=(CELLMAX-CELLMIN)/NCELL
+ CALL DOCELLS(NCELL,DCELL,CELLMIN,XCELL,YCELL)
+ ENDIF
+C
+ ANGLEI=(PI/180.)*PAR(IPPAR1(11))
+ COSI=COS(ANGLEI)
+ SINI=SIN(ANGLEI)
+ IF(ABS(COSI).LT.R4PREC) COSI=0.0
+ IF(ABS(SINI).LT.R4PREC) SINI=0.0
+C
+ CALL LUCE1(NPRINT,PRINT6,PRINTFILE,NFILE,NPRINTL,COSI,SINI,
+ 1 R4PREC,CELLMAX,CELLMIN,DCELL,NCELL,XCELL,YCELL,ICELL,
+ 2 AUS,AUS1,AUS2,AUS3,IFLAG(4))
+C
+C Normalizing to unit max. the light curve,
+C for bolometric flux and each used color band
+ K1=IPPAR1(12)+2
+ ANORMC=PAR(K1)
+ IF(ANORMC.GT.0.0) CALL NORMC(ANORMC,IFLAG(4))
+C
+C Light to observer elapsed time computation
+ IF(IFLAG(20).GT.0) THEN
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPLUC=TEMPLUC+TEMPD
+ IF(PRINT6) WRITE(6,8005) TEMPD,TEMPLUC
+C IF(PRINTFILE) WRITE(NFILE,8005)TEMPD,TEMPLUC
+ 8005 FORMAT(/' Elapsed time for received ligth :',G20.5,5X,
+ 1 ' Total:',G20.5/)
+ ENDIF
+C
+C ------------
+C =============== => | PRINTING |
+C ------------
+ IF (IFLAG(20).GT.0) THEN
+C
+ CALL PRINTING(AMOUNT,NFILE6,NFILE,IPPAR1,IPPAR2,
+ 1 NAMFLAG,DESCFLAG,NAMPAR,DESCPAR,AUS)
+C
+C I/O Elapsed time computation
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPIO=TEMPIO+TEMPD
+ IF(PRINT6) WRITE(6,8001) TEMPD,TEMPIO
+C IF(PRINTFILE) WRITE(NFILE,8001)TEMPD,TEMPIO
+C
+ ENDIF
+C
+C
+C ------------
+C =============== => | PLOTTING |
+C ------------
+ IF (IFLAG(21).GT.0) THEN
+C
+ CALL PLOTTING(PRINTFILE,NFILE,PRINT6,
+ 1 AMOUNTP,NFILE6P,NFILEP,NASCIIP)
+C
+C I/O Elapsed time computation
+ TEMPO=SECNDS(TEMPO0)
+ TEMPD=TEMPO-TEMPREC
+ TEMPREC=TEMPO
+ TEMPIOP=TEMPIOP+TEMPD
+ IF(PRINT6) WRITE(6,8006) TEMPD,TEMPIOP
+C IF(PRINTFILE) WRITE(NFILE,8006)TEMPD,TEMPIOP
+ 8006 FORMAT(/' Elapsed time for plot:',12X,G20.5,5X,' Total:',G20.5/)
+C
+ ENDIF
+C
+C
+ IFLAG(15)=-1
+ IF(MODFLAG(1).GT.0) IFLAG(1)=-1
+ IF(MODFLAG(2).GT.0) IFLAG(2)=-1
+ IF(MODFLAG(3).GT.0) IFLAG(3)=-1
+ IF(MODFLAG(12).GT.0) IFLAG(12)=-1
+ IF(MODFLAG(9).GT.0) IFLAG(9)=-1
+ IF(MODFLAG(10).GT.0) IFLAG(10)=-1
+C
+C ...........................................Return to input
+ GOTO 12
+ END
+C
+ SUBROUTINE ALLSEL(K3,NALL1,NALL2)
+C ---------------------------------------------------------
+C This routine, called by subroutine reflect, check the
+C alignement list, from NALL1 to NALL2 ,deleting all alignement
+C passing throught the object K3.
+C ---------------------------------------------------------
+ PARAMETER (MAXPT=200000 ,MAXALLIN=250000, MAXBANDE=5)
+ PARAMETER (MAXPTC=50000)
+C
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+ COMMON /ALLIN/ NALLMX,NALL,ITERDONE,
+ 1 ISOUR(MAXALLIN),JRIC(MAXALLIN),
+ 2 TRANSFI(MAXALLIN),TRANSFJ(MAXALLIN),
+ 3 FSOURI(MAXALLIN),FSOURJ(MAXALLIN),
+ 4 FRICI(MAXALLIN),FRICJ(MAXALLIN),
+ 5 RIJ(MAXALLIN),
+ 6 RXIJ(MAXALLIN),RYIJ(MAXALLIN),RZIJ(MAXALLIN),
+ 7 COSGI(MAXALLIN),COSGJ(MAXALLIN),
+ 8 IAUS(MAXALLIN)
+C
+C ..........................................................
+ LOGICAL IAUS,CONDX,CONDY,CONDZ
+C
+C to avoid looping with no points (nifc(1,k3)=nifc(2,k3)=0)
+ IF(NP123C(K3).LE.0) RETURN
+ IF(NALL1.GT.NALL2) RETURN
+C
+ DO 5 I=NALL1,NALL2
+ 5 IAUS(I)=.FALSE.
+C
+C Loops on alignements, looking for poits in K3 belonging
+C to the alignement line.
+ N1=NIFC(1,K3)
+ N2=NIFC(2,K3)
+C
+ DO 10 I=NALL1,NALL2
+ DO 20 K=N1,N2
+C
+C ........ First test into or outside interval
+C Si puo' semplificare questo groviglio di if ?
+ XIK=( XC(K) - XC(ISOUR(I)) )
+C If X outside interval Rxij
+ ABSD=ABS(DIMCX(K))
+ IF(RXIJ(I).GE.0.) THEN
+ IF(XIK .GT. RXIJ(I)+ABSD .OR. XIK .LT. -ABSD)
+ 1 GOTO 20
+ ELSE
+ IF(XIK .LT. RXIJ(I)-ABSD .OR. XIK .GT. ABSD)
+ 1 GOTO 20
+ ENDIF
+C If Y outside interval Ryij
+ YIK=( YC(K) - YC(ISOUR(I)) )
+ ABSD=ABS(DIMCY(K))
+ IF(RYIJ(I).GE.0.) THEN
+ IF(YIK .GT. RYIJ(I)+ABSD .OR. YIK .LT. -ABSD)
+ 1 GOTO 20
+ ELSE
+ IF(YIK .LT. RYIJ(I)-ABSD .OR. YIK .GT. ABSD)
+ 1 GOTO 20
+ ENDIF
+C If Z outside interval Rzij
+ ZIK=( ZC(K) - ZC(ISOUR(I)) )
+ ABSD=ABS(DIMCZ(K))
+ IF(RZIJ(I).GE.0.0) THEN
+ IF(ZIK .GT. RZIJ(I)+ABSD .OR. ZIK .LT. -ABSD)
+ 1 GOTO 20
+ ELSE
+ IF(ZIK .LT. RZIJ(I)-ABSD .OR. ZIK .GT. ABSD)
+ 1 GOTO 20
+ ENDIF
+C
+C In these cases the x,y, or z condition hasn't been tested before
+ IF(ABS(RXIJ(I)).GT.0.0) THEN
+ CONDX=.TRUE.
+ ELSE
+ CONDX=.FALSE.
+ ENDIF
+ IF(ABS(RYIJ(I)).GT.0.0) THEN
+ CONDY=.TRUE.
+ ELSE
+ CONDY=.FALSE.
+ ENDIF
+ IF(ABS(RZIJ(I)).GT.0.0) THEN
+ CONDZ=.TRUE.
+ ELSE
+ CONDZ=.FALSE.
+ ENDIF
+C ........ Then test if not on same line
+ IF(CONDX) THEN
+ IF(CONDY) THEN
+ DUM =DIMCX(K)*RYIJ(I)
+ DUM1=DIMCY(K)*RXIJ(I)
+ AMX=MAX(ABS(DUM+DUM1),ABS(DUM-DUM1))
+ IF( ABS( RYIJ(I)*XIK-RXIJ(I)*YIK ) .GT. AMX ) GOTO 20
+ ENDIF
+ IF(CONDZ) THEN
+ DUM =DIMCX(K)*RZIJ(I)
+ DUM1=DIMCZ(K)*RXIJ(I)
+ AMX=MAX(ABS(DUM+DUM1),ABS(DUM-DUM1))
+ IF( ABS( RZIJ(I)*XIK-RXIJ(I)*ZIK ) .GT. AMX ) GOTO 20
+ ENDIF
+ ENDIF
+ IF(CONDY.AND.CONDZ) THEN
+ DUM =DIMCY(K)*RZIJ(I)
+ DUM1=DIMCZ(K)*RYIJ(I)
+ AMX=MAX(ABS(DUM+DUM1),ABS(DUM-DUM1))
+ IF( ABS( RZIJ(I)*YIK-RYIJ(I)*ZIK ) .GT. AMX ) GOTO 20
+ ENDIF
+C
+C ........ The remaining surf. el. are aligned
+ IAUS(I)=.TRUE.
+C go to examine the next alignement:
+ GOTO 10
+ 20 CONTINUE
+ 10 CONTINUE
+C
+C Now alignements intercepted by K3 are deleted
+ K=NALL1-1
+ DO 30 I=NALL1,NALL2
+C If true skips over the alignement
+ IF(IAUS(I)) GOTO 30
+ K=K+1
+ IF(I.EQ.K) GOTO 30
+C If same position just increases K,otherwise move I to K
+ ISOUR(K)=ISOUR(I)
+ JRIC(K)=JRIC(I)
+ RIJ(K)=RIJ(I)
+C RXIJ,RYIJ,RZIJ are not mouved, they aren't used anymore
+ COSGI(K)=COSGI(I)
+ COSGJ(K)=COSGJ(I)
+ 30 CONTINUE
+ NALL=K
+ RETURN
+ END
+C
+ FUNCTION AMAX2(N,A)
+C ----------------------------------------------------
+C AMAX2=max(a(n))
+C ----------------------------------------------------
+ DIMENSION A(N)
+ AMAX2=A(1)
+ DO 10 I=2,N
+ IF(A(I).GT.AMAX2) AMAX2=A(I)
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE ASK(N1,N2,NAME,DESC,VALUE)
+C ----------------------------------------------------
+C Shows parameter's values and descriptors
+C ----------------------------------------------------
+ CHARACTER*(*) NAME(N2),DESC(N2)
+ DIMENSION VALUE(N2)
+ 1000 FORMAT(/(1X,A10,5X,A20,' = ',G15.7))
+ WRITE(6,1000)(NAME(J),DESC(J),VALUE(J),J=N1,N2)
+ RETURN
+ END
+C
+ SUBROUTINE ASKI(N1,N2,NAME,DESC,IVALUE)
+C ----------------------------------------------------
+C Shows flag's values and descriptors
+C ----------------------------------------------------
+ CHARACTER*(*) NAME(N2),DESC(N2)
+ DIMENSION IVALUE(N2)
+ 1000 FORMAT(/(1X,A10,5X,A20,' = ',I2))
+ WRITE(6,1000)(NAME(J),DESC(J),IVALUE(J),J=N1,N2)
+ RETURN
+ END
+C
+ SUBROUTINE BANDAC(BANDL,TEMPER,NL,DELTAL,ALAMBDA,TRASMISS)
+C ----------------------------------------------------------------
+C Computes the flux for photometric band , by integration
+C over the lambda range of (Plank function * filter transmission),
+C transmission is defined in common/filtri/.
+C Computed flux is: erg/cm**2/sec/sterad ,
+C lambda in Anstrong in /filtri/ , and in the flux computation.
+C Simpson rule is used for integration
+C BANDL= computed integral
+C alam1,alam2= LAMBDA LIMITS (UNUSED )
+C TEMPER=temperature
+C ALAMBDA(NL),TRASMISS(NL)= lambda points and their transmission coeff.
+C DELTAL=ALAMBDA(2)-ALAMBDA(1) ( Anstrong)
+C
+C ----------------------------------------------------------------
+C
+ DIMENSION ALAMBDA(NL),TRASMISS(NL)
+C
+C come e' fatta sotto e' piu' veloce
+C BANDL=0.0
+C PLAN3=PLANK(ALAMBDA(1),TEMPER)*TRASMISS(1)
+C DO 20 J=1,NL-2,2
+C PLAN1=PLAN3
+C PLAN2=PLANK(ALAMBDA(J+1),TEMPER)*TRASMISS(J+1)
+C PLAN3=PLANK(ALAMBDA(J+2),TEMPER)*TRASMISS(J+2)
+C BANDL=BANDL+(1./3.)*PLAN1+(4./3.)*PLAN2+(1./3.)*PLAN3
+C
+C DUM=+(1./3.)*PLAN1+(4./3.)*PLAN2+(1./3.)*PLAN3
+C BANDL=BANDL+DUM
+C BANDL=BANDL*DELTAL * 1.E-8
+C
+ BANDL=(1./3.)*PLANK(ALAMBDA(1),TEMPER)*TRASMISS(1)
+ DO 20 J=2,NL-3,2
+ BANDL=BANDL+(4./3.)*PLANK(ALAMBDA(J),TEMPER)*TRASMISS(J)+
+ 1 (2./3.)*PLANK(ALAMBDA(J+1),TEMPER)*TRASMISS(J+1)
+ 20 CONTINUE
+ BANDL=DELTAL * 1.E-8 *
+ 1 (BANDL+(4./3.)*PLANK(ALAMBDA(NL-1),TEMPER)*TRASMISS(NL-1)+
+ 2 (1./3.)*PLANK(ALAMBDA(NL),TEMPER)*TRASMISS(NL) )
+ RETURN
+ END
+C
+ SUBROUTINE BANDAC1(TOTTOT,TEMPER,NL,DELTAL,ALAMBDA,TRASMISS)
+C ----------------------------------------------------------------
+C UNUSED ROUTINE, THIS SHOULD SUBSTITUTE BANDAC, BEING MORE
+C PRECISE, BUT THE LINK BETWEEN PLANK LUMINOSITIES IN DIFFERENT
+C BANDS AND MAGNITUDE IS NOT WELL ESTABLISHED, FOR THIS REASON
+C YOU HAVE ALWAYS AN UNDEFINED MAG. ZERO POINT IN DIFFERENT COLORS
+C AND A DETAILED INTEGRATION AVER PLANK FUNCTION IS UNUSEFUL
+C
+C Computes the flux for photometric band , by integration
+C over the lambda range of (Plank function * filter transmission),
+C transmission is defined in common/filtri/. for few lambda points.
+C Simpson rule is used for integration, using many points and
+C interpolating (linear) over the given transmission values
+C Computed flux is: erg/cm**2/sec/sterad ,
+C lambda in Anstrong in /filtri/ , but used in cm in the
+C flux computation.
+C TOTTOT= computed integral
+C TEMPER=temperature
+C ALAMBDA(NL),TRASMISS(NL)= lambda points and their transmission coeff.
+C DELTAL = delta lambda for tabulated values in /filtri/ (UNUSED )
+C ----------------------------------------------------------------
+C
+ PARAMETER NSTEPS=20
+ DIMENSION ALAMBDA(NL),TRASMISS(NL)
+C
+ TOTTOT=0.0
+ DO 10 I=1,NL-1
+ AL1=ALAMBDA(I)
+ AL2=ALAMBDA(I+1)
+ T1=TRASMISS(I)
+ T2=TRASMISS(I+1)
+ DL=(AL2-AL1)/NSTEPS
+C changing DL from Anstrong to cm
+ DL=DL*1.E-8
+C
+ P1=PLANK(AL1,TEMPER)
+ P2=PLANK(AL2,TEMPER)
+ TOT=(1./3.)*T1*P1 + (1./3.)*T2*P2
+ DO 20 J=2,NSTEPS-2,2
+ AL=AL1+(J-1)*DL
+ T=( (T2-T1)/(AL2-AL1) ) * (AL-AL1) + T1
+ TOT=TOT + (4./3.)*T* PLANK(AL,TEMPER)
+ AL=AL1+J*DL
+ T=( (T2-T1)/(AL2-AL1) ) * (AL-AL1) + T1
+ TOT=TOT + (2./3.)* PLANK(AL,TEMPER)*T
+ 20 CONTINUE
+ AL=AL2-DL
+ T=( (T2-T1)/(AL2-AL1) ) * (AL-AL1) + T1
+ TOT=TOT + (4./3.)* PLANK(AL,TEMPER)*T
+ TOTTOT=TOTTOT+TOT*DL
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE BANDAT(NPNT,T,F,FB,FR,NL,DELTAL,ALAMBDA,TRASMISS)
+C -----------------------------------------------------------
+C This routine integrates the black body spectrum over
+C a given lambda interval for each surface element.
+C -----------------------------------------------------------
+ DIMENSION T(NPNT),F(NPNT),FB(NPNT),FR(NPNT)
+ DIMENSION ALAMBDA(NL),TRASMISS(NL)
+ PARAMETER AC4PI=1.8065E-5
+C
+ DO 10 I=1,NPNT
+ TEMPER=T(I)
+C integration over colour band
+ CALL BANDAC(TOTBAND,TEMPER,NL,DELTAL,ALAMBDA,TRASMISS)
+ F(I)=TOTBAND/(AC4PI*TEMPER**4) * (FB(I)+FR(I))
+ 10 CONTINUE
+C
+ RETURN
+ END
+C
+ SUBROUTINE BANDATM(NPNT,ALAM,T,F,FB,FR)
+C -----------------------------------------------------------
+C This routine fills F with bolometric flux FB * flux fraction
+C in a given lambda
+C FUNCTION PLANKL=PLANK*1.E-8cm = plank*delta lambda
+C -----------------------------------------------------------
+ DIMENSION T(NPNT),F(NPNT),FB(NPNT),FR(NPNT)
+ PARAMETER AC4PI=1.8065E-5
+C
+ DO 10 I=1,NPNT
+ F(I)=PLANKL(ALAM,T(I))/(AC4PI*T(I)**4) * (FB(I)+FR(I))
+ 10 CONTINUE
+C
+ RETURN
+ END
+C
+ SUBROUTINE BANDFIL(NPNT,B,FB,FR,FBANDA)
+C -----------------------------------------------------------
+C This routine assign a given fraction B of the bolometric
+C flux B to the color band FBANDA.
+C -----------------------------------------------------------
+ DIMENSION FB(NPNT),FBANDA(NPNT),FR(NPNT)
+ DO 10 I=1,NPNT
+ 10 FBANDA(I)=(FB(I)+FR(I))*B
+ RETURN
+ END
+C
+ SUBROUTINE CALCALL(ALBEDO)
+C ----------------------------------------------------
+C This routine, called by Reflect, computes the
+C remaining data for each allignement, and, for each
+C allignement I-J, inserts the inverse allignement J-I
+C ----------------------------------------------------
+ PARAMETER (MAXPT=200000 ,MAXALLIN=250000, MAXBANDE=5)
+ PARAMETER (MAXPTC=50000)
+ PARAMETER (PI4=1./(3.1415926*4.))
+ PARAMETER PI=3.1415926
+C
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+C
+ COMMON /ALLIN/ NALLMX,NALL,ITERDONE,
+ 1 ISOUR(MAXALLIN),JRIC(MAXALLIN),
+ 2 TRANSFI(MAXALLIN),TRANSFJ(MAXALLIN),
+ 3 FSOURI(MAXALLIN),FSOURJ(MAXALLIN),
+ 4 FRICI(MAXALLIN),FRICJ(MAXALLIN),
+ 5 RIJ(MAXALLIN),
+ 6 RXIJ(MAXALLIN),RYIJ(MAXALLIN),RZIJ(MAXALLIN),
+ 7 COSGI(MAXALLIN),COSGJ(MAXALLIN),
+ 8 IAUS(MAXALLIN)
+C
+C ...............................................................
+C Computes limb darkening term, allignement transfer
+C function and and set the initial reflection source
+C for each allignement
+ DO 30 K=1,NALL
+C DD= SQRT(RIJ**2)
+ DD=SQRT(RIJ(K))
+C abs(<surf.norm,rij>) -> abs(<surf.norm,unit vector rij>)=abs(cos)
+ COSGI(K)= COSGI(K) /DD
+ COSGJ(K)=-COSGJ(K) /DD
+C source and receiving surf. el. for reflection
+ KI=ISOUR(K)
+ KJ=JRIC(K)
+C albedo for source and receiving surf.el.
+ XCOSGI= 1.0 - ( XLIMBC(KI) * (1.0-COSGI(K)) )
+ XCOSGJ= 1.0 - ( XLIMBC(KJ) * (1.0-COSGJ(K)) )
+ ALBI=ALBC(KI)
+ ALBJ=ALBC(KJ)
+C transfer function from I to J surf.el.
+C TRANSFI(K)=(ALBEDO*PI4)*ALBJ*XCOSGI *COSGJ(K)*
+C 1 AC(KJ)/RIJ(K)
+C TRANSFI(K)=(ALBEDO)*ALBJ*XCOSGI *(COSGJ(K)*COSGI(K))*
+C 1 AC(KJ)/RIJ(K)/PI
+ TRANSFI(K)=(ALBEDO)*ALBJ*XCOSGI *(COSGJ(K)*COSGI(K))*
+ 1 AC(KJ)/RIJ(K)/ (PI*(1.-XLIMBC(KJ)/3.) )
+C inverse transfer from J to I
+C TRANSFJ(K)=(ALBEDO*PI4)*ALBI*XCOSGJ *COSGI(K)*
+C 1 AC(KI)/RIJ(K)
+C TRANSFJ(K)=(ALBEDO)*ALBI*XCOSGJ *(COSGJ(K)*COSGI(K))*
+C 1 AC(KI)/RIJ(K)/PI
+ TRANSFJ(K)=(ALBEDO)*ALBI*XCOSGJ *(COSGJ(K)*COSGI(K))*
+ 1 AC(KI)/RIJ(K)/ (PI*(1.-XLIMBC(KI)/3.) )
+C initial source for reflection from I to J and inverse path
+ FSOURI(K)=FBC(KI)
+ FSOURJ(K)=FBC(KJ)
+ 30 CONTINUE
+C
+ RETURN
+ END
+C
+ SUBROUTINE CANCEL(K,NFILE)
+C -----------------------------------------------------
+C In the surface element list in common /surf/
+C all the data of object K are deleted, by reordering
+C of the list.
+C ----------------------------------------------------
+ PARAMETER (MAXPT=200000 , MAXBANDE=5)
+ PARAMETER (MAXPTC=50000)
+ PARAMETER (NVECT=21+MAXBANDE)
+ PARAMETER (NVECTC=15)
+ COMMON /SURF/ NPNTMX,NPNT,NP(3),NLIM(2,3),
+ 1 VECT(MAXPT,NVECT)
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPC(3),NLIMC(2,3),
+ 1 VECTC(MAXPTC,NVECTC)
+C COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+C A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+C 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+C 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+C 5 AC(MAXPTC),
+C 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+C 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+C 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+C 8 NUMFINI(MAXPTC)
+C DIMENSION NP123C(3)
+C DIMENSION NIFC(2,3)
+C EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+C Deleting object K in common /surf/
+C
+ N1=NLIM(1,K)
+ N2=NLIM(2,K)+1
+ NPNTK=N2-N1
+ IF(NPNTK.NE.NP(K)) THEN
+ WRITE(6,1000) K,NP,NLIM
+ 1000 FORMAT(' ERROR! The fine grid surface element list descriptors',
+ 1 ' are inconsistent for object:',I5/' descriptors are:'/
+ 2 ' number of points for each object:',3I5/
+ 3 ' first and last element for each object:',6I5)
+ IF(NFILE.GT.0.AND.NFILE.LE.99) WRITE(NFILE,1000) K,NP,NLIM
+ ENDIF
+C
+ IF(N2.LE.NPNT) THEN
+C otherwise a simple change of dimension descriptors deletes the object
+ DO 10 I=1,NVECT
+ DO 10 J=N2,NPNT
+ VECT(J-NPNTK,I)=VECT(J,I)
+ 10 CONTINUE
+ DO 20 I=1,3
+ IF(I.EQ.K.OR.NLIM(1,I).LT.N2) GOTO 20
+ NLIM(1,I)=NLIM(1,I)-NPNTK
+ NLIM(2,I)=NLIM(2,I)-NPNTK
+ 20 CONTINUE
+ ENDIF
+ NLIM(1,K)=0
+ NLIM(2,K)=0
+ NP(K)=0
+ NPNT=NPNT-NPNTK
+C
+C Deleting object K in common /surfc/
+C
+C Here NLIMC(2,k) is not used, it can be zero if the
+C coarse grid map has been computed, but not the coarse surf. elem.
+C also in this case this common must be compressed; its space
+C is allocated to object K in any case.
+C
+ N1=NLIMC(1,K)
+ N2=N1+NPC(K)
+ NPNTK=N2-N1
+ IF(N2.LT.N1.OR.N2.GT.NPNTC+1.OR.
+ 1 (NLIMC(2,K).GT.0.AND.N2.NE.NLIMC(2,K)+1) ) THEN
+ WRITE(6,3000) K,NPC,NLIMC
+ 3000 FORMAT(' ERROR! The coarse grid surface element list descripto',
+ 1 'rs are inconsistent for object:',I5/' descriptors are:'/
+ 2 ' number of points for each object:',3I5/
+ 3 ' first and last element for each object:',6I5)
+ IF(NFILE.GT.0.AND.NFILE.LE.99) WRITE(NFILE,3000) K,NPC,NLIMC
+ ENDIF
+C
+ IF(N2.LE.NPNTC) THEN
+C otherwise a simple change of dimension descriptors deletes the object
+ DO 30 I=1,NVECTC
+ DO 30 J=N2,NPNTC
+ VECTC(J-NPNTK,I)=VECTC(J,I)
+ 30 CONTINUE
+ DO 40 I=1,3
+ IF(I.EQ.K.OR.NLIMC(1,I).LT.N2) GOTO 40
+ NLIMC(1,I)=NLIMC(1,I)-NPNTK
+ IF(NLIMC(2,I).GT.0) NLIMC(2,I)=NLIMC(2,I)-NPNTK
+ 40 CONTINUE
+ ENDIF
+ NLIMC(1,K)=0
+ NLIMC(2,K)=0
+ NPC(K)=0
+ NPNTC=NPNTC-NPNTK
+ RETURN
+ END
+C
+ SUBROUTINE COARSE(PRINT6,PRINTFILE,NFILE,IS)
+C ---------------------------------------------
+C From the coarse map: NUMCOARSE, containig, for
+C each fine surface element, the number of its
+C coarse surf. elem., builts the coarse surf. elem.
+C ------------------------------------------------
+C
+ PARAMETER (MAXPT=200000 , MAXBANDE=5 )
+ PARAMETER (MAXPTC=50000 )
+C
+ COMMON /SURF/ NPNTMX,NPNT,NPNT1,NPNT2,NPNT3,
+ A N1I,N1F,N2I,N2F,N3I,N3F,
+ 1 X(MAXPT),Y(MAXPT),Z(MAXPT),
+ 2 G(MAXPT),GX(MAXPT),GY(MAXPT),GZ(MAXPT),
+ 3 FX(MAXPT),FY(MAXPT),FZ(MAXPT),
+ 4 TX(MAXPT),TY(MAXPT),TZ(MAXPT),
+ 5 T(MAXPT),A(MAXPT),
+ 6 FB(MAXPT),F(MAXPT,MAXBANDE),FRIFL(MAXPT),
+ 7 XLIMB(MAXPT),NUMOBJ(MAXPT),ALB(MAXPT),
+ 8 NUMCOARSE(MAXPT)
+ DIMENSION NP123(3)
+ DIMENSION NIF(2,3)
+ EQUIVALENCE (NP123(1),NPNT1),(NIF(1,1),N1I)
+C
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+ LOGICAL PRINT6,PRINTFILE
+C
+ N1=NIF(1,IS)
+ N2=NIF(2,IS)
+C
+ N1C=NIFC(1,IS)
+ N2C=N1C+NP123C(IS)-1
+C Nifc(2,.)=0 when coarse map NUMCOARSE exist, but NOT surf.el.
+ NIFC(2,IS)=N2C
+C
+ DO 10 I=N1C,N2C
+ XC(I)=0.0
+ YC(I)=0.0
+ ZC(I)=0.0
+ GXC(I)=0.0
+ GYC(I)=0.0
+ GZC(I)=0.0
+ AC(I)=0.0
+ FBC(I)=0.0
+ FRIFLC(I)=0.0
+ XLIMBC(I)=0.0
+ ALBC(I)=0.0
+ DIMCX(I)=0.0
+ DIMCY(I)=0.0
+ DIMCZ(I)=0.0
+ NUMFINI(I)=0
+ NUMOBJC(I)=IS
+ 10 CONTINUE
+C
+ DO 20 I=N1,N2
+ K=NUMCOARSE(I)
+ IF(K.GT.N2C.OR.K.LT.N1C) THEN
+ WRITE(6,1000) I,NUMOBJ(I),K,IS
+ IF(PRINTFILE) WRITE(NFILE,1000) I,NUMOBJ(I),K,IS
+ 1000 FORMAT(/' ERROR! , surface element:',I6,' belonging to the'
+ 1 ' object:',I3/' assigned to coarse point:',I6,
+ 2 ' belonging to the wrong object:',I3)
+ GOTO 20
+ ENDIF
+C The coarse X,Y,Z etc are obtained by a weighted mean of
+C their fine surf.interval (fine surf.el. areas are weights);
+C to shorten the calculation you can simply sum The fine's x,y,z;
+C then divide by NUMFINI. This could be a not too bad approximation.
+ NUMFINI(K)=NUMFINI(K)+1
+ XC(K)=XC(K)+X(I) *A(I)
+ YC(K)=YC(K)+Y(I) *A(I)
+ ZC(K)=ZC(K)+Z(I) *A(I)
+ GXC(K)=GXC(K)+GX(I) *A(I)
+ GYC(K)=GYC(K)+GY(I) *A(I)
+ GZC(K)=GZC(K)+GZ(I) *A(I)
+ AC(K)=AC(K)+A(I)
+ FBC(K)=FBC(K)+FB(I)
+ XLIMBC(K)=XLIMBC(K)+XLIMB(I) *A(I)
+ ALBC(K)=ALBC(K)+ALB(I) *A(I)
+ DIMCX(K)=DIMCX(K)+ FX(I)+TX(I)
+ DIMCY(K)=DIMCY(K)+ FY(I)+TY(I)
+ DIMCZ(K)=DIMCZ(K)+ FZ(I)+TZ(I)
+ 20 CONTINUE
+C
+ DO 30 I=N1C,N2C
+ XC(I)=XC(I) /AC(I)
+ YC(I)=YC(I) /AC(I)
+ ZC(I)=ZC(I) /AC(I)
+ GXC(I)=GXC(I)/AC(I)
+ GYC(I)=GYC(I)/AC(I)
+ GZC(I)=GZC(I)/AC(I)
+ XLIMBC(I)=XLIMBC(I)/AC(I)
+ ALBC(I)=ALBC(I) /AC(I)
+ 30 CONTINUE
+C renormalizing G
+ DO 35 I=N1C,N2C
+ GC(I)=SQRT(GXC(I)*GXC(I)+GYC(I)*GYC(I)+GZC(I)*GZC(I))
+ GXC(I)=GXC(I)/GC(I)
+ GYC(I)=GYC(I)/GC(I)
+ GZC(I)=GZC(I)/GC(I)
+ 35 CONTINUE
+C
+ RETURN
+ END
+C
+ SUBROUTINE CYL(PRINT6,PRINTFILE,NFILE,NRF,NTHF,NCORD,PKRI,NOBJ,
+ 1 NPUNT,NPUNTC,R,RI,H,HI,X0,Y0,Z0,ALZO,ROTAZ,R4PREC,
+ 2 NPNTMX,X,Y,Z,G,GX,GY,GZ,FX,FY,FZ,TX,TY,TZ,AREA,A,NUMCOARSE,
+ 3 NUMOBJ, R2,RR,RFACTOR,SINF,COSF)
+C ----------------------------------------------------------------
+C Drawn a cylindrica surface, see sfera2 for general comments.
+C NRF= number of fine r mesh in which R-RI is divided.
+C NTHF= coarsing factor=number of fine for each coarse
+C NRC= number of r coarse
+C NFC= number of phi coarse (r dependent and >=4 )
+C After drawing on the plane x-y the cylinder can be rotated:
+C ROTAZ= anti-clockwise angle of rotazion on the plane (x-x'angle)
+C ALZO= clockwise rotation around the y' axis (z-z" angle)
+C Ater the rotation the disk is shifted from (0,0,0) in (x0,y0,z0)
+C
+C NOBJ: number of this object, numobj(side surf.el)=nobj or nobj+3
+C if angle>r4prec ( to represent a concave obj. by two finctious ones,
+C avoiding in this way disk transparency in luce1 )
+C ######## se si abbandona la simmetria di rotazione per il disco
+C si devono aumentare il numero di oggetti che lo rappresentano, perche'
+C non diventi trasparente ########################################
+C ----------------------------------------------------------------
+ DIMENSION X(NPNTMX),Y(NPNTMX),Z(NPNTMX)
+ DIMENSION G(NPNTMX),GX(NPNTMX),GY(NPNTMX),GZ(NPNTMX)
+ DIMENSION FX(NPNTMX),FY(NPNTMX),FZ(NPNTMX)
+ DIMENSION TX(NPNTMX),TY(NPNTMX),TZ(NPNTMX)
+ DIMENSION A(NPNTMX)
+ DIMENSION NUMCOARSE(NPNTMX),NUMOBJ(NPNTMX)
+ DIMENSION R2(NTHF),RR(NTHF),RFACTOR(NTHF),SINF(NTHF),COSF(NTHF)
+ LOGICAL PRINTFILE,PRINT6
+ DATA PI /3.1415926/
+C DATA PI /3.141592653589793/
+C REAL*8 PI
+C
+ WIDTH5=(H-HI)*0.5
+ HI2=HI*0.5
+C Disk thikening angle
+ ANGLE=ATAN2(WIDTH5,R-RI)
+ COSA=COS(ANGLE)
+ SINA=SIN(ANGLE)
+ IF(ABS(SINA).LT.R4PREC) SINA=0.0
+ TANA=TAN(ANGLE)
+ IF(ABS(TANA).LT.R4PREC) TANA=0.0
+ COSA1=1./COSA
+ NOBJ1=NOBJ
+ IF(ANGLE.GT.R4PREC) NOBJ1=NOBJ+3
+C
+C Coarse mesh
+ NRC=NRF/NTHF
+ IF(MOD(NRF,NTHF).NE.0) THEN
+ WRITE(6,1000) NRF,NTHF
+ 1000 FORMAT(' ERROR! In drawing a cylinder the number of r mesh:',
+ 1 I5,' is not consistent with the coarsing factor:',I5)
+ ENDIF
+C
+C 1/cosa isn't included here: meshing is done on the disk central plane
+ DRC=(R-RI)/NRC
+ DRC5=DRC*0.5
+C fine r mesh
+ DRF=DRC/NTHF
+ DRF5=DRF*0.5
+ DRF5TANA=DRF5*TANA
+C
+C NPUNTC will be incremented below. (NPUNTC is the first void coarse)
+C IC is instead the actual coarse, or the last filled one.
+C NPUNT is the number of surf.el. generated by this routine
+ IC=NPUNTC-2
+ NPUNT=0
+ AREA=0.0
+C
+C ................... Coarse R loop ( from center RI to R)
+ DO 10 IRC=1,NRC
+ R1C=RI+(IRC-1)*DRC
+ R2C=R1C+DRC
+C Number of phi coarse for this Radius
+ NFC=(2.*PI)*R2C/DRC
+ IF(NFC.LT.4) NFC=4
+ DFC=(2.*PI)/NFC
+C
+C Radius fine for this R coarse
+ DUM=R1C
+ DO 15 I=1,NTHF
+ R1 =DUM
+ R2(I)=DUM+DRF
+ RR(I)=DUM+DRF5
+ RFACTOR(I)=COSA1 * (R2(I)-R1)*(R2(I)+R1)
+ DUM=DUM+DRF
+ 15 CONTINUE
+C
+C ...................... PHI coarse loop
+ DO 10 IFC=1,NFC
+ F1C=(IFC-1)*DFC
+ F2C=F1C+DFC
+C coarse surf.el. number
+ IC=IC+2
+C PHI fine interval
+ DFF=DFC/NTHF
+ DFF5=DFF*0.5
+C Cord, arch or tangent value for surf el. extension
+ IF(NCORD.EQ.2) THEN
+ SINF5=SIN(DFF5)
+ ELSE IF (NCORD.GE.3) THEN
+ TANF5=TAN(DFF5)
+ ENDIF
+C
+C First phi fine value fi coarse refers to the interval boundaries,
+C fi fine is shifted at the interval center (+DFF5)
+ F=F1C-DFF5
+C sin and cos(f) values for this range if phi fine
+ DO 35 I=1,NTHF
+ F=F+DFF
+ COSF(I)=COS(F)
+ SINF(I)=SIN(F)
+ 35 CONTINUE
+C
+C .................... LOOP on phi fine interv. of coarse
+ DO 10 IFF=1,NTHF
+C
+C .................... LOOP on R fine interval of coarse
+ DO 10 IRF=1,NTHF
+C
+ IF(NCORD.LE.1) THEN
+ CORDA=DFF5 * R2(IRF)
+ ELSE IF(NCORD.EQ.2) THEN
+ CORDA=SINF5 * R2(IRF)
+ ELSE
+ CORDA=TANF5 * R2(IRF)
+ ENDIF
+C
+C Upper surface point
+ IF(NPUNT.GE.NPNTMX) THEN
+ LOST=LOST+2
+ GOTO 10
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)= RR(IRF)*COSF(IFF)
+ Y(NPUNT)= RR(IRF)*SINF(IFF)
+ Z(NPUNT)= HI2 + RR(IRF)*TANA
+ GX(NPUNT)=-SINA*COSF(IFF)
+ GY(NPUNT)=-SINA*SINF(IFF)
+ GZ(NPUNT)= COSA
+ FX(NPUNT)=-CORDA*SINF(IFF)
+ FY(NPUNT)= CORDA*COSF(IFF)
+ FZ(NPUNT)=0.
+ TX(NPUNT)= DRF5*COSF(IFF)
+ TY(NPUNT)= DRF5*SINF(IFF)
+ TZ(NPUNT)= DRF5TANA
+ A(NPUNT)=DFF5*RFACTOR(IRF)
+ NUMCOARSE(NPUNT)=IC
+ NUMOBJ(NPUNT)=NOBJ
+C area computation
+ AREA=A(NPUNT)+AREA
+C
+C Lower surface point
+ NPUNT1=NPUNT
+ NPUNT=NPUNT+1
+ X(NPUNT)= X(NPUNT1)
+ Y(NPUNT)= Y(NPUNT1)
+ Z(NPUNT)=-Z(NPUNT1)
+ GX(NPUNT)= GX(NPUNT1)
+ GY(NPUNT)= GY(NPUNT1)
+ GZ(NPUNT)=-GZ(NPUNT1)
+ FX(NPUNT)=FX(NPUNT1)
+ FY(NPUNT)=FY(NPUNT1)
+ FZ(NPUNT)=0.
+ TX(NPUNT)= TX(NPUNT1)
+ TY(NPUNT)= TY(NPUNT1)
+ TZ(NPUNT)= TZ(NPUNT1)
+ A(NPUNT)=A(NPUNT1)
+ NUMCOARSE(NPUNT)=IC+1
+ NUMOBJ(NPUNT)=NOBJ
+C Total area computation
+ AREA=A(NPUNT)+AREA
+C
+ 10 CONTINUE
+C
+C The IC+1 coarse has been filled by loop 10
+ IC=IC+1
+C
+C If this cylinder had a widht the side surface is drawn
+C
+ IF(H.LE.R4PREC) GOTO 400
+C
+C coarse and fine h and fi interval definition
+ NHC=INT(H/DRC)
+ IF(NHC.LE.0) NHC=1
+ DHC=H/NHC
+ DHC5=DHC*0.5
+ DHF=DHC/NTHF
+ DHF5=DHF*0.5
+ NFC=(2.*PI)*R/DRC
+ IF(NFC.LT.4) NFC=4
+ DFC=(2.*PI)/NFC
+ DFF=DFC/NTHF
+ DFF5=DFF*0.5
+C
+ IF(NCORD.LE.1) THEN
+ CORDA=R*DFF5
+ ELSE IF(NCORD.EQ.2) THEN
+ CORDA=R*SIN(DFF5)
+ ELSE
+ CORDA=R*TAN(DFF5)
+ ENDIF
+C
+ AREAR=R*DFF*DHF
+C
+C ........................ LOOP on phi coarse
+ DO 60 IFC=1,NFC
+ F1C=(IFC-1)*DFC
+C sinf and cosf fine of this coarse
+ DUM=F1C-DFF5
+ DO 65 I=1,NTHF
+ DUM=DUM+DFF
+ COSF(I)=COS(DUM)
+ SINF(I)=SIN(DUM)
+ 65 CONTINUE
+C
+C ....................... Loop on H coarse
+ DO 60 IHC=1,NHC
+C lower z value (going from -H to (h-dhc))
+ H1=-H*0.5-DHF5+(IHC-1)*DHC
+C
+ IC=IC+1
+C
+C ...................... LOOP on phi fine
+ DO 60 IFF=1,NTHF
+C
+C ........................ LOOP on H fine
+ DO 60 IHF=1,NTHF
+ HF=H1 + IHF*DHF
+C
+ IF(NPUNT.GE.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 60
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)=R*COSF(IFF)
+ Y(NPUNT)=R*SINF(IFF)
+ Z(NPUNT)=HF
+ GX(NPUNT)=COSF(IFF)
+ GY(NPUNT)=SINF(IFF)
+ GZ(NPUNT)=0.0
+ FX(NPUNT)=-CORDA*SINF(IFF)
+ FY(NPUNT)= CORDA*COSF(IFF)
+ FZ(NPUNT)=0.
+ TX(NPUNT)=0.0
+ TY(NPUNT)=0.0
+ TZ(NPUNT)= -DHF5
+ A(NPUNT)=AREAR
+ NUMCOARSE(NPUNT)=IC
+ NUMOBJ(NPUNT)=NOBJ1
+C Total area computation
+ AREA=AREA + A(NPUNT)
+C
+ 60 CONTINUE
+C
+ 400 CONTINUE
+C internal surface is drawn
+ IF(HI.LE.R4PREC.OR.RI.LE.R4PREC.OR.PKRI.LE.0) GOTO 500
+C
+C coarse and fine height mesh
+ NHC=INT(HI/DRC)
+ IF(NHC.LE.0) NHC=1
+ DHC=HI/NHC
+ DHC5=DHC*0.5
+ DHF=DHC/NTHF
+ DHF5=DHF*0.5
+C phi coarse mesh
+ NFC=(2.*PI)*RI/DRC
+ IF(NFC.LT.4) NFC=4
+ DFC=(2.*PI)/NFC
+ DFF=DFC/NTHF
+ DFF5=DFF*0.5
+C Area for each fine interval
+ AREAR=RI*DFF*DHF
+C
+ IF(NCORD.LE.1) THEN
+ CORDA=RI*DFF5
+ ELSE IF(NCORD.EQ.2) THEN
+ CORDA=RI*SIN(DFF5)
+ ELSE
+ CORDA=RI*TAN(DFF5)
+ ENDIF
+C ........................ LOOP on phi coarse
+ DO 70 IFC=1,NFC
+ F1C=(IFC-1)*DFC
+C sinf and cosf fine of this coarse
+ DUM=F1C-DFF5
+ DO 75 I=1,NTHF
+ DUM=DUM+DFF
+ COSF(I)=COS(DUM)
+ SINF(I)=SIN(DUM)
+ 75 CONTINUE
+C
+C ....................... Loop on H coarse
+ DO 70 IHC=1,NHC
+C lower z value (going from -Hi to (hi-dhc))
+ H1=-HI*0.5-DHF5 + (IHC-1)*DHC
+C
+ IC=IC+1
+C
+C ...................... LOOP on phi fine
+ DO 70 IFF=1,NTHF
+C
+C ........................ LOOP on H fine
+ DO 70 IHF=1,NTHF
+ HF=H1+ IHF*DHF
+C
+ IF(NPUNT.GE.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 70
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)=RI*COSF(IFF)
+ Y(NPUNT)=RI*SINF(IFF)
+ Z(NPUNT)=HF
+ GX(NPUNT)=COSF(IFF)
+ GY(NPUNT)=SINF(IFF)
+ GZ(NPUNT)=0.0
+ FX(NPUNT)=-CORDA*SINF(IFF)
+ FY(NPUNT)= CORDA*COSF(IFF)
+ FZ(NPUNT)=0.
+ TX(NPUNT)=0.0
+ TY(NPUNT)=0.0
+ TZ(NPUNT)=DHF5
+ A(NPUNT)=AREAR
+ NUMCOARSE(NPUNT)=IC
+ NUMOBJ(NPUNT)=NOBJ1
+C Total area computation
+ AREA=A(NPUNT)+AREA
+C
+ 70 CONTINUE
+C .................. END OF DISK SURFACE DRAWING !
+ 500 CONTINUE
+ NPUNTC=IC+1
+C
+C Constant potential USEFUL FOR DISK ?
+ DO 59 I=1,NPUNT
+ 59 G(I)=1.
+C
+ IF(LOST.GT.0) THEN
+ WRITE(6,2000) NPUNT,LOST
+ IF(PRINTFILE) WRITE(NFILE,1000) NPUNT,LOST
+ 2000 FORMAT(' ERROR! In Cyl:OBJECT NOT PROPERLY DRAWN! '/
+ 1 ' space available for only:',I6,' surface points.',
+ 2 ' lost points:',I6/
+ 3 ' REDUCE THE NPHI PARAMETER OR GO INTO THE FORTRAN LIST TO',
+ 4 ' INCREASE MAXPT')
+ ENDIF
+C
+C ................ optional disk rotation .....
+C
+ IF(ALZO.EQ.0.0.AND.ROTAZ.EQ.0.0) GOTO 800
+C The disk is rotated by rotaz around the z axis (anti-clockwise)
+C then around the new Y axis (clockwise) by alzo
+ COSP=COS(ROTAZ)
+ SINP=SIN(ROTAZ)
+ COSI=COS(ALZO)
+ SINI=SIN(ALZO)
+C
+ IF(ABS(COSP).LT.R4PREC) COSP=0.0
+ IF(ABS(SINP).LT.R4PREC) SINP=0.0
+ IF(ABS(COSI).LT.R4PREC) COSI=0.0
+ IF(ABS(SINI).LT.R4PREC) SINI=0.0
+C
+ DO 80 I=1,NPUNT
+ X1= X(I)* (COSP*COSI) + Y(I) * SINP - Z(I) * (COSP*SINI)
+ Y1=-X(I)* (SINP*COSI) + Y(I) * COSP + Z(I) * (SINP*SINI)
+ Z1= X(I)* SINI + Z(I) * COSI
+ X(I)=X1
+ Y(I)=Y1
+ Z(I)=Z1
+ X1= GX(I)* (COSP*COSI) + GY(I) * SINP - GZ(I) * (COSP*SINI)
+ Y1=-GX(I)* (SINP*COSI) + GY(I) * COSP + GZ(I) * (SINP*SINI)
+ Z1= GX(I)* SINI + GZ(I) * COSI
+ GX(I)=X1
+ GY(I)=Y1
+ GZ(I)=Z1
+ X1= TX(I)* (COSP*COSI) + TY(I) * SINP - TZ(I) * (COSP*SINI)
+ Y1=-TX(I)* (SINP*COSI) + TY(I) * COSP + TZ(I) * (SINP*SINI)
+ Z1= TX(I)* SINI + TZ(I) * COSI
+ TX(I)=X1
+ TY(I)=Y1
+ TZ(I)=Z1
+ X1= FX(I)* (COSP*COSI) + FY(I) * SINP - FZ(I) * (COSP*SINI)
+ Y1=-FX(I)* (SINP*COSI) + FY(I) * COSP + FZ(I) * (SINP*SINI)
+ Z1= FX(I)* SINI + FZ(I) * COSI
+ FX(I)=X1
+ FY(I)=Y1
+ FZ(I)=Z1
+C
+ 80 CONTINUE
+C
+ 800 CONTINUE
+C shifting the disk from the origin=(0,0,0) to x0,y0,z0
+ IF(X0.NE.0.0) THEN
+ DO 90 I=1,NPUNT
+ 90 X(I)=X(I)+X0
+ ENDIF
+ IF(Y0.NE.0.0) THEN
+ DO 91 I=1,NPUNT
+ 91 Y(I)=Y(I)+Y0
+ ENDIF
+ IF(Z0.NE.0.0) THEN
+ DO 92 I=1,NPUNT
+ 92 Z(I)=Z(I)+Z0
+ ENDIF
+C
+ RETURN
+ END
+C
+ SUBROUTINE DOCELLS(N,D,DMIN,X,Y)
+C -----------------------------------------
+C Computes the grid cell CENTERS for the
+C projection plane
+C ----------------------------------------
+ DIMENSION X(N),Y(N)
+ DO 10 I=1,N
+ X(I)=(DMIN+D*0.5)+(I-1)*D
+ Y(I)=-X(I)
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE FASGRID(N,FASE)
+C -----------------------------------------------------
+C Defines an equispaced set of phases at which the
+C light to observer will be computed
+C -----------------------------------------------------
+ DIMENSION FASE(N)
+ REAL*8 DEL
+ DEL=(6.283185307179586)/N
+ DO 10 I=1,N
+ 10 FASE(I)=DEL*(I-1)
+ RETURN
+ END
+C
+ SUBROUTINE FASREAD(NFASI,NFILEF,NFILE,MAX,FASE)
+C ---------------------------------------------------
+C Reads fase values from unit NFILEF
+C ---------------------------------------------------
+ DIMENSION FASE(MAX)
+ DATA PI2/6.2831853/
+ DO 10 I=1,MAX
+ READ(NFILEF,*,END=500) FASE(I)
+ IF(FASE(I).LT.0.0.OR.FASE(I).GT.1.0) GOTO 500
+ 10 CONTINUE
+ WRITE(6,1000)MAX,MAX
+ IF(NFILE.GT.0.AND.NFILE.LE.99)WRITE(NFILE,1000)MAX,MAX
+ 1000 FORMAT(' WARNING! You have reached the maximum number of',
+ 1 ' phases: ',I6/' If you need a grater value go into',
+ 2 ' the fortran list to change MAXFASI=',I6)
+ 500 CONTINUE
+ NFASI=I-1
+ IF(NFASI.LE.0) THEN
+ WRITE(6,1001) NFILEF
+ IF(NFILE.GT.0.AND.NFILE.LE.99)WRITE(NFILE,1001)NFILEF
+ 1001 FORMAT(' ERROR! Zero phase values read from unit:',I5)
+ RETURN
+ ENDIF
+ DO 20 I=1,NFASI
+ 20 FASE(I)=FASE(I)*PI2
+ RETURN
+ END
+C
+ SUBROUTINE FILFILT
+C ------------------------------------------------------
+C filling common /FILTRI/ with lambda and transmission
+C factors for Johnson bands UBVRI .
+C ( data from C.W.Allen - Astrophysical Quantities -
+C Athone Press - London 1976 - pg.201 )
+C - - - - - - - - - - - - - - - - - -
+C Lambda values must be equi-spaced
+C - - - - - - - - - - - - - - - - - -
+C Number of lambda values must be odd
+C - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+C PARAMETERS :=NU,NB,NV,NR,NI = num lambda for each filter,(<=MAXLAM)
+C NLAMBDA(maxfilt)=NU,NB,NV,NR,NI = num lambda for each filter,
+C NLAMBDA must be <= nfiltl
+C DELTALAM= delta lambda for each filter
+C U(nu),B(nb),..... are the light fraction tranmitted for lambdas
+C in UL(nu),BL(nb) .... (in Anstrong ) ,
+C LAMBDA in common /filtri/ are Anstrong,
+C ------------------------------------------------------
+C
+ PARAMETER (NU=9 , NB=11 , NV=13, NR=1, NI=1 )
+C
+ PARAMETER ( MAXFILT=5 , MAXLAM=15 )
+ COMMON /FILTRI/ NFILT,NLAM, NLAMBDA(MAXFILT),DELTALAM(MAXLAM),
+ 1 ALAMBDA(MAXLAM,MAXFILT),TRASMISS(MAXLAM,MAXFILT),
+ 2 COMPFRAC(MAXLAM,3)
+C
+ DIMENSION U(NU),B(NB),V(NV),R(NR),AI(NI)
+ DIMENSION UL(NU),BL(NB),VL(NV),RL(NR),AIL(NI)
+C
+ DATA U/0.0,0.13,0.6,0.92,1.,0.72,0.07,0.0,0.0/
+ DATA UL/2800.,3000.,3200.,3400.,3600.,3800.,4000.,4200.,4400./
+ DATA B/0.0,0.13,0.92,1.00,0.92,0.76,0.56,0.39,0.2,0.07,0.0/
+ DATA BL/3600.,3800.,4000.,4200.,4400.,4600.,4800.,5000.,
+ 1 5200.,5400.,5600./
+ DATA V/0.0,0.01,0.36,0.91,0.98,0.8,0.59,0.39,0.22,0.09,
+ 1 0.03,0.01,0.0/
+ DATA VL/4600.,4800.,5000.,5200.,5400.,5600.,5800.,6000.,6200.,
+ 1 6400.,6600.,6800.,7000./
+C
+C NFILT=MAXFILT
+C NLAM=MAXLAM
+C
+C Number of lambda for each band
+C
+ NLAMBDA(1)=NU
+ NLAMBDA(2)=NB
+ NLAMBDA(3)=NV
+ NLAMBDA(4)=NR
+ NLAMBDA(5)=NI
+C
+ DELTALAM(1)=200.
+ DELTALAM(2)=200.
+ DELTALAM(3)=200.
+ DELTALAM(4)=200.
+ DELTALAM(5)=200.
+C
+C .............................................................
+C Filling of common /filtri/
+C The following is not dependent on lambda and transmission
+C .............................................................
+C
+C following loop to avoid bad changes to this routine,
+C by inexpert users inserting their preferred photometric bands.
+ DO 5 I=1,NFILT
+ IF(NLAMBDA(I).GT.NLAM) THEN
+ WRITE(6,1000) NLAMBDA(I),I,NLAM
+ 1000 FORMAT(' ERROR ! : WRONG CHANGE IN SUBROUTINE FILFILT'/
+ 1 1X,I5,'=NUMBER OF LAMBDA FOR BAND',I2,' > MAXLAM=',I5)
+ NLAMBDA(I)=NLAM
+ ENDIF
+ 5 CONTINUE
+C
+C U - BAND - Filter 1
+C
+ DO 10 I=1,NU
+ ALAMBDA(I,1)=UL(I)
+ 10 TRASMISS(I,1)=U(I)
+C
+C
+C B - BAND - Filter 2
+C
+ DO 20 I=1,NB
+ ALAMBDA(I,2)=BL(I)
+ 20 TRASMISS(I,2)=B(I)
+C
+C
+C V - BAND - Filter 3
+C
+ DO 30 I=1,NV
+ ALAMBDA(I,3)=VL(I)
+ 30 TRASMISS(I,3)=V(I)
+C
+C
+C R - BAND - Filter 4
+C
+ DO 40 I=1,NR
+ ALAMBDA(I,4)=RL(I)
+ 40 TRASMISS(I,4)=R(I)
+C
+C
+C I - BAND - Filter 5
+C
+ DO 50 I=1,NI
+ ALAMBDA(I,5)=AIL(I)
+ 50 TRASMISS(I,5)=AI(I)
+C
+ RETURN
+ END
+C
+ SUBROUTINE FILL(N,A,B)
+C -----------------------------
+C Fills B(i) with A
+C ----------------------------
+ DIMENSION B(N)
+ DO 10 I=1,N
+ 10 B(I)=A
+ RETURN
+ END
+C
+ SUBROUTINE FILL1(N,A,B,C)
+C ----------------------------
+C Fills B(i) with A*C(i)
+C ----------------------------
+ DIMENSION B(N),C(N)
+ DO 10 I=1,N
+ 10 B(I)=C(I)*A
+ RETURN
+ END
+C
+ SUBROUTINE FINDALL(K1,K2,R4PREC)
+C ----------------------------------------------
+C Finds allignement between the objects K1 and K2.
+C An allignement is a straight line between a
+C surface element in k1 and one in k2 which
+C doesn't intersect the objects k1 or k2.
+C This way of finding alignements doesn't work
+C for complicated surfaces like ones with holes
+C and mountains, but only when the outward normal
+C line to a surface element never meets the surface:
+C here a path between the two object is
+C recognized as an alignement when the two surface
+C elements at its ends look each other.
+C
+C Sin and cos < r4prec are assumed 0 ( sin,cos(0) is not precise)
+C r4prec is inserted to avoiding double precision here.
+C ----------------------------------------------
+ PARAMETER (MAXPT=200000 ,MAXALLIN=250000, MAXBANDE=5)
+ PARAMETER (MAXPTC=50000)
+C
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+C
+ COMMON /ALLIN/ NALLMX,NALL,ITERDONE,
+ 1 ISOUR(MAXALLIN),JRIC(MAXALLIN),
+ 2 TRANSFI(MAXALLIN),TRANSFJ(MAXALLIN),
+ 3 FSOURI(MAXALLIN),FSOURJ(MAXALLIN),
+ 4 FRICI(MAXALLIN),FRICJ(MAXALLIN),
+ 5 RIJ(MAXALLIN),
+ 6 RXIJ(MAXALLIN),RYIJ(MAXALLIN),RZIJ(MAXALLIN),
+ 7 COSGI(MAXALLIN),COSGJ(MAXALLIN),
+ 8 IAUS(MAXALLIN)
+C
+C ................................................................
+ N1=NIFC(1,K1)
+ N2=NIFC(2,K1)
+ N3=NIFC(1,K2)
+ N4=NIFC(2,K2)
+ LOST=0
+C looks for alignements
+C ----------------------------------- LOOP BEGINS
+ DO 10 I=N1,N2
+ DO 10 J=N3,N4
+C Vector joining the two surface elements (from I to J)
+ RIJX=XC(J)-XC(I)
+ RIJY=YC(J)-YC(I)
+ RIJZ=ZC(J)-ZC(I)
+C surface normal in I projected along this vector
+ COSGI1=GXC(I)*RIJX+GYC(I)*RIJY+GZC(I)*RIJZ
+C If < 0 the surface element I doesn't look toward J
+ IF (COSGI1.LE.R4PREC) GOTO 10
+C IF (COSGI1.LE.0.) GOTO 10
+C surface normal in J projected along vector from I to J
+ COSGJ1=GXC(J)*RIJX+GYC(J)*RIJY+GZC(J)*RIJZ
+C If > 0 the surface element J doesn't look toward I
+ IF (COSGJ1.GE.-R4PREC) GOTO 10
+C IF (COSGJ1.GE.0.) GOTO 10
+C Here the path is recognized as an alignement
+C Test if the alignement list is full.
+ IF(NALL.GT.NALLMX) THEN
+ LOST=LOST+1
+ GOTO 10
+ ENDIF
+C Insert this path in the alignement list
+ NALL=NALL+1
+ ISOUR(NALL)=I
+ JRIC(NALL)=J
+ RIJ(NALL)=RIJX*RIJX+RIJY*RIJY+RIJZ*RIJZ
+ RXIJ(NALL)=RIJX
+ RYIJ(NALL)=RIJY
+ RZIJ(NALL)=RIJZ
+ COSGI(NALL)=COSGI1
+ COSGJ(NALL)=COSGJ1
+ 10 CONTINUE
+C ------------------------- LOOP ENDS
+ IF(LOST.GT.0) THEN
+ WRITE(6,1000) LOST
+ 1000 FORMAT(' ERROR: ',I6,' ALIGNEMENTS LOST IN REFLECTION '/
+ 1 ' Reduce the number of surface elements or go into the'/
+ 2 ' FORTRAN list to change parameter MAXALLIN')
+ ENDIF
+ RETURN
+ END
+C
+ SUBROUTINE FINDCOMM(N1,N2,NAME,COMMAND,KF)
+C -----------------------------------------------------
+C Looks for "COMMAND" in a list of NAMEs,
+C Returns KF = position of COMMAND in vector NAME
+C returns KF=0 if COMMAND isn't found in NAME vactor
+C -----------------------------------------------------
+ CHARACTER*(*) NAME(N2),COMMAND
+ DO 10 I=N1,N2
+ IF(NAME(I).NE.COMMAND) GOTO 10
+ KF=I
+ RETURN
+ 10 CONTINUE
+ KF=0
+ RETURN
+ END
+C
+ SUBROUTINE FINEFRIFL
+C ------------------------------------------------
+C From coarse grid surface elements reflection source
+C computes fine grid surf. elem. rreflect. source
+C ---------------------------------------------------
+C
+ PARAMETER (MAXPT=200000 , MAXBANDE=5 )
+ PARAMETER (MAXPTC=50000 )
+C
+ COMMON /SURF/ NPNTMX,NPNT,NPNT1,NPNT2,NPNT3,
+ A N1I,N1F,N2I,N2F,N3I,N3F,
+ 1 X(MAXPT),Y(MAXPT),Z(MAXPT),
+ 2 G(MAXPT),GX(MAXPT),GY(MAXPT),GZ(MAXPT),
+ 3 FX(MAXPT),FY(MAXPT),FZ(MAXPT),
+ 4 TX(MAXPT),TY(MAXPT),TZ(MAXPT),
+ 5 T(MAXPT),A(MAXPT),
+ 6 FB(MAXPT),F(MAXPT,MAXBANDE),FRIFL(MAXPT),
+ 7 XLIMB(MAXPT),NUMOBJ(MAXPT),ALB(MAXPT),
+ 8 NUMCOARSE(MAXPT)
+ DIMENSION NP123(3)
+ DIMENSION NIF(2,3)
+ EQUIVALENCE (NP123(1),NPNT1),(NIF(1,1),N1I)
+C
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+ DO 10 I=1,NPNT
+ K=NUMCOARSE(I)
+ IF(FRIFLC(K).LE.0.0) THEN
+ FRIFL(I)=0.0
+ ELSE
+ FRIFL(I)=A(I)/AC(K)*FRIFLC(K)
+ ENDIF
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE GRADOMEGA(X,Y,Z,GX,GY,GZ,G,Q)
+C ---------------------------------------------------
+C Used by subr. Roches : normalized omega gradient
+C ---------------------------------------------------
+ SQ =1./(X*X+Y*Y+Z*Z)**1.5
+ SQ1=1./((1.-X)*(1.-X)+Y*Y+Z*Z)**1.5
+ A= SQ+SQ1*Q
+ B= A -1.-Q
+ GX= X * B + Q * (1.-SQ1)
+ GY= Y * B
+ GZ= Z * A
+ G=SQRT(GX*GX+GY*GY+GZ*GZ)
+ GX=GX/G
+ GY=GY/G
+ GZ=GZ/G
+ RETURN
+ END
+C
+ SUBROUTINE GRAVITY1(KK,NPNT,BETA,TINPUT,G,T)
+C -------------------------------------------------
+C Compute the star surface temperature variation
+C due to the different gravity. (gravity darkening)
+C BETA is the gravity darkening exponent
+C TINPUT the max surface temperature of the star
+C as given in input.(= pole temperature)
+C -------------------------------------------------
+ DIMENSION G(NPNT),T(NPNT)
+C
+C Temperature profile computation
+C and max temperature of the computed profile
+ T(1)=1.
+ TMAX=1.
+ KK=1
+ DO 10 I=2,NPNT
+ T(I)=(G(I)/G(1))**BETA
+ IF(T(I).GT.TMAX) THEN
+ TMAX=T(I)
+ KK=I
+ ENDIF
+ 10 CONTINUE
+C Renormalizes the temp. profile to give the input medium temp.
+ DO 20 I=1,NPNT
+ T(I)=(TINPUT/TMAX)*T(I)
+ 20 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE INPUT(NAMFLAG,DESCFLAG,IFLAGDEF,NAMPAR,DESCPAR,
+ 1 PARDEF,IPPAR1,IPPAR2,STOP)
+C ------------------------------------------------------------
+C Reads input and sets flags and parameters
+C ------------------------------------------------------------
+ PARAMETER (MAXFLAG=21 , MAXPAR=128 , MAXPARS=20 )
+ PARAMETER (MAXTITLE=5)
+ DIMENSION IPPAR1(MAXFLAG),IPPAR2(MAXFLAG),MODFLAG(MAXFLAG)
+ DIMENSION IFLAGDEF(MAXFLAG),PARDEF(MAXPAR)
+ CHARACTER *10 NAMFLAG(MAXFLAG)
+ CHARACTER *20 DESCFLAG(MAXFLAG)
+ CHARACTER *10 NAMPAR(MAXPAR)
+ CHARACTER *20 DESCPAR(MAXPAR)
+ LOGICAL STOP
+C
+ CHARACTER*80 COMMAND,COMMAND1
+ CHARACTER PROMPT *15, PROMPT1 *20
+ DATA PROMPT,PROMPT1/' INPUT FLAG >',' INPUT PARAMETER >'/
+C
+ COMMON /PAR/ NFLAG,IFLAG(MAXFLAG),NPAR,PAR(MAXPAR)
+C DIMENSION PARS(MAXPARS,3)
+C EQUIVALENCE (PARS(1,1),PAR(1))
+ COMMON/TITLE/NTITLEMX,NTITLE,TITLE(MAXTITLE)
+ CHARACTER*80 TITLE
+C
+C ------------------------------------------------ INPUT FLAG LOOP
+ GOTO 12
+C ........................... Loop entry for input: "?"
+ 8 CALL ASKI(1,NFLAG,NAMFLAG,DESCFLAG,IFLAG)
+ GOTO 12
+C ........................... Loop entry for bad command
+ 10 WRITE(6,1000) COMMAND
+1000 FORMAT(' FLAG NAME NOT RECOGNIZED: ',A10,' Please reenter')
+C ........................... Loop entry for normal reading
+ 12 CALL UPPERC(COMMAND,PROMPT,KFLU)
+C If input "?" loops again
+ IF(KFLU.NE.0) GOTO 8
+C Looks for flag name
+ CALL FINDCOMM(1,NFLAG,NAMFLAG,COMMAND,KFL)
+C If name not found loops again
+ IF(KFL.LE.0.OR.KFL.GT.NFLAG) GOTO 10
+C If iflag(19)=1 set flag off
+ IF(IFLAG(19).EQ.1) THEN
+ IFLAG(KFL)=0
+ WRITE(6,1100) KFL,NAMFLAG(KFL),DESCFLAG(KFL)
+ 1100 FORMAT(' Set off: flag number:',I5,5x,A10,2X,A20)
+ GOTO 12
+ ENDIF
+C Otherwise set flag on
+C (16 is the entry from parameter loop)
+ 16 IFLAG(KFL)=1
+ WRITE(6,1110) KFL,NAMFLAG(KFL),DESCFLAG(KFL)
+ 1110 FORMAT(' Set on: flag number:',I5,5x,A10,2X,A20)
+C Exits from input loops
+ IF(IFLAG(15).EQ.1) GOTO 300
+ IF(IFLAG(16).EQ.1.OR.IFLAG(17).EQ.1) GOTO 500
+C If there aren't flag's parameters loops again
+ IF(IPPAR1(KFL).LE.0) GOTO 12
+C
+C ------------------------------------------ INPUT PARAMETERS LOOP
+ GOTO 22
+C ........................... Loop entry for input: "?"
+ 18 CALL ASK(IPPAR1(KFL),IPPAR2(KFL),NAMPAR,DESCPAR,PAR)
+ GOTO 22
+C ........................... Loop entry for bad command
+ 20 WRITE(6,2000) COMMAND1
+2000 FORMAT(' PARAMETER NOT RECOGNIZED: ',A10,' Please reenter')
+C ........................... Loop entry for normal reading
+ 22 CALL UPPERC(COMMAND1,PROMPT1,KFLU1)
+C If input "?" loops again
+ IF(KFLU1.NE.0) GOTO 18
+C Looks for parameter name
+ CALL FINDCOMM(IPPAR1(KFL),IPPAR2(KFL),NAMPAR,COMMAND1,KFL1)
+C If name is found sets its value
+ IF(KFL1.GE.IPPAR1(KFL).AND.KFL1.LE.IPPAR2(KFL)) THEN
+ GOTO 24
+ 23 WRITE(6,2001)
+ 2001 FORMAT(' INPUT ERROR! REENTER')
+ 24 WRITE(6,2002)
+C Following format is not ANSI standard
+ 2002 FORMAT(' VALUE > ',$)
+ READ(5,*,ERR=23,END=23) COMM
+ 2003 FORMAT(E20.0)
+ PAR(KFL1)=COMM
+ WRITE(6,2100) KFL1,NAMPAR(KFL1),DESCPAR(KFL1),PAR(KFL1)
+ 2100 FORMAT(' Parameter number:',I5,5x,A10,2X,A20,' = ',G15.7)
+ ELSE
+C If not found tests if it's a flag
+ CALL FINDCOMM(1,NFLAG,NAMFLAG,COMMAND1,KFL2)
+C Not being a flag loops again
+ IF(KFL2.LE.0.OR.KFL2.GT.NFLAG) GOTO 20
+C Otherwise goes to flag loop
+ KFL=KFL2
+ GOTO 16
+ ENDIF
+ GOTO 22
+C ------------------------------------------ END OF PARAMETER LOOP
+C ------------------------------------------ END OF FLAG LOOP
+C
+ 300 STOP=.FALSE.
+ RETURN
+ 500 STOP=.TRUE.
+ RETURN
+ END
+C
+ SUBROUTINE INULL(N,IA)
+C -----------------------
+C Fills IA with 0 .
+C Depending on computers,
+C NULL could be also used.
+C ----------------------
+ DIMENSION IA(N)
+ DO 10 I=1,N
+ 10 IA(I)=0
+ RETURN
+ END
+C
+ SUBROUTINE LAGR(Q,PRECISION)
+C ---------------------------------------------------------
+C Lagrange point computation for the mass ratio Q
+C the distance between the two body is the measure unit(=1)
+C PRECISION = accurancy of Newton-Rapson algorithm,
+C remember that you need 1.E-7 to have all the 7 digits of a real*4 number
+C ---------------------------------------------------------
+ COMMON /ROCHE/ AL1,AL2,AL3,XA1,XA2,XB1,XB2,YA,YB,
+ 1 OMEGAL1,OMEGAL2,OMEGAL3,
+ 2 RPOLE(3),XPOLE(3),YPOLE(3),ZPOLE(3),TPOLE(3),
+ 3 GXPOLE(3),GYPOLE(3),GZPOLE(3)
+C
+ EXTERNAL LAGRF1,LAGRF2,LAGRF3
+ OM(Q,XV)=1./ABS(XV)+Q/ABS(1.-XV)+(1.+Q)/2.*XV*XV-Q*XV
+C
+C Initial values ( used if the computation below doesn't succeed)
+C
+ AL1=0.0
+ AL2=0.0
+ AL3=0.0
+C Approx transv. dimension for critical lobes(Plavec BAC,15:165(64))
+C Plavec formula assume q<1 :
+ IF(Q.LE.1) THEN
+ ALOGQ=LOG10(Q)
+ ELSE
+ ALOGQ=LOG10(1./Q)
+ ENDIF
+C
+ YA=10**(-0.4221-0.2084*ALOGQ)
+ YB=10**(-0.4236+0.2743*ALOGQ)
+C Lagrange points: Inner point L1
+C L2 (near minor body) R1,R2 are search limits for rtsafe routine
+C L3 (near major body )
+ NERROR=0
+ R1=0.0
+ R2=1.0
+ AL1=RTSAFE(LAGRF1,R1,R2,PRECISION,Q,DUM,NERROR)
+ OMEGAL1=OM(Q,AL1)
+ R1=1.
+ R2=2.
+ AL2=RTSAFE(LAGRF2,R1,R2,PRECISION,Q,DUM,NERROR)
+ OMEGAL2=OM(Q,AL2)
+ R1=-1.
+ R2=0.
+ AL3=RTSAFE(LAGRF3,R1,R2,PRECISION,Q,DUM,NERROR)
+ OMEGAL3=OM(Q,AL3)
+C
+ IF(NERROR.GT.0) THEN
+ WRITE(6,1000)
+ 1000 FORMAT(/' ERROR In subr. LAGR, Incorrect Lagrange Points')
+ ENDIF
+ RETURN
+ END
+C
+ FUNCTION LAGRF1(X,F,DF,Q,DUM)
+C -----------------------------------------------------------------
+C Intersection x-axis <-> roche lobes for a given Q, called by LAGR
+C near L1 point
+C -----------------------------------------------------------------
+C F= X**5*(1.+Q)- X**4*(2.+3.*Q)+X**3*(1.+3.*Q)-X**2+X*2.-1.
+C DF=5.*X**4*(1.+Q)-4.*X**3*(2.+3.*Q)+3.*X**2*(1.+3.*Q)-2.*X +2.
+ Q1=1.+Q
+ QXQ=X*Q1-Q
+ QXQ2=2.*QXQ
+ RX=1.-X
+ RX2=RX*RX
+ X2=X*X
+ F=RX2 *( X2 * QXQ -1. ) + X2*Q
+ DF= RX* ( RX* X* (QXQ2 + X*Q1) - X * (X*QXQ2) -2. ) +2.*X*Q
+ RETURN
+ END
+C
+ FUNCTION LAGRF2(X,F,DF,Q,DUM)
+C -----------------------------------------------------------------
+C Intersection x-axis <-> roche lobes for a given Q, called by LAGR
+C near L2 point
+C -----------------------------------------------------------------
+ Q1=1.+Q
+ QXQ=X*Q1-Q
+ QXQ2=2.*QXQ
+ RX=X-1.
+ RX2=RX*RX
+ X2=X*X
+ F=RX2 *( X2 * QXQ -1. ) - X2*Q
+ DF= RX* ( RX* X* (QXQ2 + X*Q1) + X * (X*QXQ2) -2. ) -2.*X*Q
+ RETURN
+ END
+C
+ FUNCTION LAGRF3(X,F,DF,Q,DUM)
+C -----------------------------------------------------------------
+C Intersection x-axis <-> roche lobes for a given Q, called by LAGR
+C near L3 point
+C -----------------------------------------------------------------
+ Q1=1.+Q
+ QXQ=X*Q1-Q
+ QXQ2=2.*QXQ
+ RX=1.-X
+ RX2=RX*RX
+ X2=X*X
+ F=RX2 *( X2 * QXQ +1. ) + X2*Q
+ DF= RX* ( RX* X* (QXQ2 + X*Q1) + X * (X*QXQ2) +2. ) +2.*X*Q
+ RETURN
+ END
+C
+ SUBROUTINE LBOLT(NPNT,F,T,A)
+C -----------------------------------------------------------
+C This routine computes the bolometric flux integrated over
+C each surface element, assuming a black body emission
+C flux=erg/cm**2/sec/sterad .=caT**4/(4 PI)
+C F=FLUX*AREA=erg/sec/sterad
+C we could set, to shorten calculation: flux=T**4, being a*c/(4*pi)=1.
+C -----------------------------------------------------------
+ DIMENSION F(NPNT),T(NPNT),A(NPNT)
+ PARAMETER ( AC4PI=7.567E-15 *3.E10/4./3.1415926 )
+ DO 10 I=1,NPNT
+ F(I)=T(I)**4*A(I) * AC4PI
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE LIMB(N,D,A,T,XLIMB,ALB)
+C -----------------------------------------------
+C D is the input given limb darkening; A the albedo
+C A time dependece could be inserted here,
+C but now this is the same as the Fill routine.
+C -----------------------------------------------
+ DIMENSION T(N),XLIMB(N),ALB(N)
+ DO 10 I=1,N
+ ALB(I)=A
+ 10 XLIMB(I)=D
+ RETURN
+ END
+C
+ SUBROUTINE LOBES(Q,PRECISION,OMEGA,R,IS)
+C -----------------------------------------------------------------
+C Given Q and a potential value (omega), computes the
+C Intersections of the x-axis with Roche lobes with potential OMEGA
+C and mass ratio Q : IS=1 for the major star, 2 for the minor one
+C If Omega is not given (<= 0.0) then it is deduced from R :
+C intersection of lobe and X axis
+C Precision= precision in Newton-Cotes routine. 1.E-7 to give all the 7
+C digits of a real*4 number
+C -----------------------------------------------------------------
+ COMMON/ROCHE/ AL1,AL2,AL3,XA1,XA2,XB1,XB2,YA,YB,
+ 1 OMEGAL1,OMEGAL2,OMEGAL3,
+ 2 RPOLE(3),XPOLE(3),YPOLE(3),ZPOLE(3),TPOLE(3),
+ 3 GXPOLE(3),GYPOLE(3),GZPOLE(3)
+ EXTERNAL LOBESF1,LOBESF2,LOBESF3
+ LOGICAL PRIMO
+C
+C OM(Q,XV)=1./XV+Q/(1.-XV)+(1.+Q)/2.*XV*XV-Q*XV+Q*Q/2./(1.+Q) old omega def.
+ OM(Q,XV)=1./ABS(XV)+Q/ABS(1.-XV)+(1.+Q)/2.*XV*XV-Q*XV
+C
+ PRIMO=.TRUE.
+ IF(IS.EQ.3) THEN
+ WRITE(6,1000)
+ 1000 FORMAT(' ERROR ! In subroutine LOBES: object 3 can''t be a',
+ 1 ' Roche model star!')
+ RETURN
+ ENDIF
+C
+ 10 IF(OMEGA.LE.0.0) THEN
+ WRITE(6,2000) R
+ 2000 FORMAT(' WARNING ! : Subr.Lobes computes Omega from R:',
+ 1 G15.7/' R is intended as the lobe radius at the ',
+ 2 ' intersection with x axis, near L1' )
+ IF(R.LE.0..OR.R.GT.1.0) THEN
+ IF(IS.EQ.1) THEN
+ WRITE(6,2010) R,AL1
+ R=AL1
+ ELSE
+ WRITE(6,2010) R,1.-AL1
+ R=1.-AL1
+ ENDIF
+ 2010 FORMAT(' WARNING ! : Incorrect R :',G14.7,
+ 1 ' Setting critical lobe: R=',G14.7)
+ ENDIF
+ IF(IS.EQ.1.) THEN
+ OMEGA=OM(Q,R)
+ ELSE
+ OMEGA=OM(Q,1.-R)
+ ENDIF
+ WRITE(6,2020) OMEGA
+ 2020 FORMAT(' I set OMEGA=',G15.7)
+ PRIMO=.FALSE.
+ ENDIF
+C R1,R2 =bounds for searching the inner lobe intersection
+C R3,R4 for outher lobe intersetion
+C For major star:
+ IF(IS.EQ.1) THEN
+ 20 R1=0.0
+ R2=AL1
+ R3=-AL1
+ R4=0.
+C Error counter :
+ NERROR=0
+ XA1=RTSAFE(LOBESF1,R1,R2,PRECISION,Q,OMEGA,NERROR)
+ XA2=RTSAFE(LOBESF3,R3,R4,PRECISION,Q,OMEGA,NERROR)
+C
+C Error recovery: 1: try to compute omega from R ( if not done)
+C 2: set R=critical lobe dimension (if not done)
+C 3: as 2 and guess values for XA1,XA2=star bounds
+ IF(NERROR.GT.0.OR.XA1.GT.AL1.OR.XA2.LT.AL3) THEN
+ WRITE(6,3000) IS
+ 3000 FORMAT(' ERROR ! wrong Q or OMEGA ; object:',I3)
+ WRITE(6,4000)AL1,AL2,AL3,XA1,XA2,OMEGA
+ 4000 FORMAT(' Lagrange points:',3G15.7/
+ 1 ' Surf.- X axis intersections:',2G14.7,' Omega:',G14.7)
+C
+ IF(PRIMO) THEN
+ IF(R.LE.0.) THEN
+ WRITE(6,2010) R,AL1
+ R=AL1
+ ENDIF
+ OMEGA=OM(Q,R)
+ WRITE(6,2000) R
+ WRITE(6,2020) OMEGA
+ PRIMO=.FALSE.
+ GOTO 20
+ ELSE IF(R.NE.AL1) THEN
+ WRITE(6,2010) R,AL1
+ R=AL1
+ OMEGA=OM(Q,R)
+ WRITE(6,2000) R
+ WRITE(6,2020) OMEGA
+ GOTO 20
+ ELSE
+ WRITE(6,5000) XA1,XA2,AL1
+ 5000 FORMAT(' ERROR RECOVERY NOT POSSIBLE:',
+ 1 ' wrong intersections:',2G15.7/
+ 2 ' used L1:',G14.7)
+ XA1=AL1
+ XA2=-AL1
+ ENDIF
+ ENDIF
+ ELSE
+C Minor star
+ 30 R1=AL1
+ R2=1.
+ R3=1.
+ R4=AL2
+C Error counter reset
+ NERROR=0
+ XB1=RTSAFE(LOBESF1,R1,R2,PRECISION,Q,OMEGA,NERROR)
+ XB2=RTSAFE(LOBESF2,R3,R4,PRECISION,Q,OMEGA,NERROR)
+C Error recovery
+ IF(NERROR.GT.0.OR.XB1.LT.AL1.OR.XB2.GT.AL2) THEN
+ WRITE(6,3000) IS
+ WRITE(6,4000)AL1,AL2,AL3,XB1,XB2,OMEGA
+C
+ IF(PRIMO) THEN
+ IF(R.LE.0.) THEN
+ WRITE(6,2010) R,1.-AL1
+ R=1.-AL1
+ ENDIF
+ OMEGA=OM(Q,1.-R)
+ WRITE(6,2000)R
+ WRITE(6,2020) OMEGA
+ PRIMO=.FALSE.
+ GOTO 30
+ ELSE IF(R.NE.1.-AL1) THEN
+ WRITE(6,2010) R,1.-AL1
+ R=1.-AL1
+ OMEGA=OM(Q,AL1)
+ WRITE(6,2000)R
+ WRITE(6,2020) OMEGA
+ GOTO 20
+ ELSE
+ WRITE(6,5000) XB1,XB2,AL1
+ XB1=AL1
+ XB2=1.+AL1
+ ENDIF
+ ENDIF
+ ENDIF
+ RETURN
+ END
+C
+ SUBROUTINE LOBESF1(X,F,DF,Q,OMEGA)
+C ------------------------------------------------
+C used by lobes : intersection of roche lobe
+C and X axis near L1 for a given omega and q
+C ------------------------------------------------
+ A=(1.+Q)/2.
+ B=A+Q
+C C=OMEGA-Q- (Q**2/2./(1.+Q)) old omega definition
+C D=Q-1.-OMEGA+(Q**2/2./(1.+Q))
+ C=OMEGA-Q
+ D=Q-1.-OMEGA
+ F=-A*X**4+B*X**3+C*X**2+D*X+1.
+ DF=-4.*A*X**3+3.*X**2*B+2.*C*X+D
+ RETURN
+ END
+C
+ SUBROUTINE LOBESF2(X,F,DF,Q,OMEGA)
+C ------------------------------------------------
+C used by lobes : intersection of roche lobe
+C and X axis near L2(after star B)for a given omega and q
+C ------------------------------------------------
+ A=(1.+Q)/2.
+ B=A+Q
+C C=OMEGA-Q- (Q**2/2./(1.+Q)) old omega definition
+C D= - Q-1.-OMEGA+(Q**2/2./(1.+Q))
+ C=OMEGA-Q
+ D= - Q-1.-OMEGA
+ F=-A*X**4+B*X**3+C*X**2+D*X+1.
+ DF=-4.*A*X**3+3.*X**2*B+2.*C*X+D
+ RETURN
+ END
+C
+ SUBROUTINE LOBESF3(X,F,DF,Q,OMEGA)
+C ------------------------------------------------
+C used by lobes : intersection of roche lobe
+C and X axis near L3(before star A)for a given omega and q
+C ------------------------------------------------
+ A=(1.+Q)/2.
+ B=A+Q
+C C=OMEGA-Q- (Q**2/2./(1.+Q)) old omega definition
+C D=Q+1.-OMEGA+(Q**2/2./(1.+Q))
+ C=OMEGA-Q
+ D=Q+1.-OMEGA
+ F=-A*X**4+B*X**3+C*X**2+D*X -1.
+ DF=-4.*A*X**3+3.*X**2*B+2.*C*X+D
+ RETURN
+ END
+C
+ SUBROUTINE LUCE1(NPRINT,PRINT6,PRINTFILE,NFILE,NPRINTL,
+ 1 COSI,SINI,R4PREC,CELLMAX,CELLMIN,DCELL,NCELL,
+ 2 XCELL,YCELL,ICELL,COSG,NEXT,ILOST,IGAINED,IBAND)
+C ------------------------------------------------------------
+C This routine computes the light received by the observer
+C at each phase. The image of the system is projected
+C onto a plane perpendicular to the observer at each phase.
+C This plane is divided into a number of grid points.
+C Two special cases are to be considered:
+C a) overfill case: more surf. elements of the same object
+C fall into the same grid cell,( this happens at the
+C star disk boundary); in this case a linklist is
+C constructed with these points.
+C b) underfill case: a surface element's projection is greater
+C than the cell dimension. To account for this the projected
+C surf.el. boundary are computed and the surf.el. is expanded
+C over its cells on the grid.
+C If two surface elements of different objects falls into the
+C same cell, the further's light is decreased by the cell contribution,
+C the part of its light falling into the cell being eclipsed.
+C
+C Finally, an integration over the perpendicular plane gives
+C the total received light for the phase.
+C
+C COS < R4PREC is assumed 0, this shorten alittle bit the computations,
+C by eliminating uneffective points, seen side-face. Points at 90 degrees
+C are made visible by rounding errors in cos(about 90) whitout this trick.
+C
+C XCELL,YCELL are the center of the cells, used only for printing
+C The projection reference system is LEFT-HANDED (for better
+C vision of objects on the screen),the other systems are right-handed!
+C YCELL is built in DOCELLS going from ymax to -ymax,
+C making apparence of the printed projection plane right-handed again.
+C #################################### NOTA #################
+C La distanza dall'osservatore e' ricalcolata ogni volta;
+C tenerla in un vettore dist risparmia di doverla ricalcolare quando
+C serve, ma ogni volta che si deve usare si deve fare un if per vedere
+C se e' gia' calcolata oppure si deve calcolare. Per questo l'uso
+C di dist va provato in vari casi, per vedere se conviene.
+C #############################################################
+C
+C PRINTFILE : if true can print information on unit NFILE on each
+C NPRINT phase steps; NPRINTL defines the amount of printing
+C PRINT6 : if true some log printing on unit 6
+C -------------------------------------------------------------
+ PARAMETER (MAXPT=200000 , MAXFASI=1000, MAXBANDE=5)
+C
+ COMMON /SURF/ NPNTMX,NPNT,NPNT1,NPNT2,NPNT3,
+ A N1I,N1F,N2I,N2F,N3I,N3F,
+ 1 X(MAXPT),Y(MAXPT),Z(MAXPT),
+ 2 G(MAXPT),GX(MAXPT),GY(MAXPT),GZ(MAXPT),
+ 3 FX(MAXPT),FY(MAXPT),FZ(MAXPT),
+ 4 TX(MAXPT),TY(MAXPT),TZ(MAXPT),
+ 5 T(MAXPT),A(MAXPT),
+ 6 FB(MAXPT),F(MAXPT,MAXBANDE),FRIFL(MAXPT),
+ 7 XLIMB(MAXPT),NUMOBJ(MAXPT),ALB(MAXPT),
+ 8 NUMCOARSE(MAXPT)
+ DIMENSION NP123(3)
+ DIMENSION NIF(2,3)
+ EQUIVALENCE (NP123(1),NPNT1),(NIF(1,1),N1I)
+C
+ COMMON /LIGHT/ NFASIMX,NFASI,NBANDE,
+ 1 FASE(MAXFASI),O(MAXBANDE,MAXFASI),
+ 2 C(3,MAXBANDE,MAXFASI),CT(MAXBANDE,MAXFASI),
+ 3 CBOL(3,MAXFASI),CBOLT(MAXFASI),
+ 4 CBOLR(3,MAXFASI),CBOLTR(MAXFASI),
+ 5 AREA(3),FBOLTOT(3),FBOL(3),FBOLREF(3)
+C
+ DIMENSION XCELL(NCELL),YCELL(NCELL),ICELL(NCELL,NCELL)
+C
+ DIMENSION IBAND(MAXBANDE)
+ DIMENSION COSG(MAXPT),NEXT(MAXPT),IGAINED(MAXPT),ILOST(MAXPT)
+C Print buffer, used to store data to be printed in the same line
+ PARAMETER (MAXBUF=20)
+ DIMENSION NBUF(MAXBUF)
+ LOGICAL PRINTFILE,PRINT6
+ DIMENSION XESP12(4)
+ EQUIVALENCE (XESP1,XESP12(1)),(XESP2,XESP12(2))
+ EQUIVALENCE (XESP3,XESP12(3)),(XESP4,XESP12(4))
+C ............................................................
+ CALL NULL(MAXBANDE*MAXFASI*3,C)
+ CALL NULL(MAXBANDE*MAXFASI,CT)
+ CALL NULL(MAXFASI*3,CBOL)
+ CALL NULL(MAXFASI,CBOLT)
+ CALL NULL(MAXFASI*3,CBOLR)
+ CALL NULL(MAXFASI,CBOLTR)
+C ------------------------------------- Loop on phases
+ IF(PRINT6) TEMPO0=SECNDS(0.0)
+ DO 10 IF=1,NFASI
+C ................. time printing for previous phase
+ IF(PRINT6) THEN
+ TEMPO2=SECNDS(TEMPO0)
+ WRITE(6,8999) IF,FASE(IF),TEMPO2
+ TEMPO0=TEMPO0+TEMPO2
+ 8999 FORMAT(' Sub.LUCE1: beginning phase:',I5,' : ',G12.5,
+ 1 ' elapsed time for last phase:',G15.7)
+ ENDIF
+C .................. counters for statistic
+C ######### poi da eliminare, sprecano tempo! ############
+ NCOUNTS=0 ! considered surf. el.
+ NCOUNTL=0 ! surf. el put in link list
+ NCOUNTE=0 ! cell expansions
+ NCOUNTV=0 ! " " over void cells
+ NCOUNTD=0 ! " " over different object's
+ NCOUNTD1=0 ! " " (of which eclipsing)
+ NCOUNTD2=0 ! " " (of which eclipsed )
+ NCOUNTF=0 ! " " over same object's
+ NCOUNTT=0 ! total remaining surf points
+ NCOUNTTL=0 ! remaining surf points in link lists
+ NCOUNTL1=0 ! number of remaining link lists
+C
+ COSF=COS(FASE(IF))
+ SINF=SIN(FASE(IF))
+C Components of the unit vector toward the observer Z0
+C and of unit vectors of the perpendicular plane X0,Y0
+C X0,Y0,Z0 = Y , -Z , X to make printing easier and nicer.
+C In this way we have a left-handed projection reference system
+C and a right-handed object reference system,
+C but, when printing the projection plane, the Y coordinate
+C is printed as -Y, making this printing right-handed again.
+ Y0X=COSI*COSF
+ Y0Y=-COSI*SINF
+ Y0Z=-SINI
+ X0X=SINF
+ X0Y=COSF
+ X0Z=0.0
+ Z0X=COSF*SINI
+ Z0Y=-SINF*SINI
+ Z0Z=COSI
+C
+C tempo=secnds(0.0)
+ CALL INULL(NCELL*NCELL,ICELL)
+C tempo1=secnds(tempo)
+C write(6,9000) if,tempo1
+C if(printfile) write(nfile,9000) if,tempo1
+ 9000 format(' Sub.luce1: phase:',I5,
+ 1 ' elapsed time for grid initialization:',G15.7)
+C
+C ------------------------------------ Loop on surface elements
+ LOST=0
+ DO 30 I=1,NPNT
+C Number of cells covered by the surf. elem.
+ IGAINED(I)=0
+C normal to the surface element projected toward the observer
+ DUM = Z0X*GX(I)+Z0Y*GY(I)+Z0Z*GZ(I)
+C Skips the element if it doesn't look toward the observer
+ IF(DUM.LE.R4PREC) GOTO 30
+C
+ NCOUNTS=NCOUNTS+1
+ COSG(I)=DUM
+ NEXT(I)=0
+C Number of cells covered by the object, but lost being eclipsed
+ ILOST(I)=0
+C
+C Looks for the point of the plane in which the element falls
+C In the following statement X0z=0. CAN BE is assumed ####
+ XX0=X(I)*X0X+Y(I)*X0Y+Z(I)*X0Z -CELLMIN
+ YY0=X(I)*Y0X+Y(I)*Y0Y+Z(I)*Y0Z -CELLMIN
+C find the cell in which the point falls
+ IX=INT(XX0/DCELL)+1
+ IY=INT(YY0/DCELL)+1
+C Set the upper boundary point (useful for xx0==cellmax)
+ IF(IX.EQ.NCELL+1) IX=NCELL
+ IF(IY.EQ.NCELL+1) IY=NCELL
+C Safety test
+ IF(IX.LE.0.OR.IX.GT.NCELL.OR.IY.LE.0.OR.IY.GT.NCELL) THEN
+ LOST=LOST+1
+ GOTO 30
+ ENDIF
+C ............... Expand each projected surf.point on the grid
+C to its projected dimension; decreasing light
+C of partially eclipsing surface elements if
+C belonging to different objects.
+C ################### nota #####################################
+C FPX,FPY,TPX,TPY,YMIN,YMAX,YESP possono tutti essere shiftati
+C di YY0, invece di partire da 0. Idem per XPX,XPY etc, si
+C risparmia tempo
+C ##############################################################
+C Computing the projected dimensions for the i surf. element:
+C F and T are projected onto the visual plane
+C ( In the following statement X0z=0. CAN BE is assumed ####)
+ FPX= FX(I)*X0X + FY(I)*X0Y + FZ(I)*X0Z
+ FPY= FX(I)*Y0X + FY(I)*Y0Y + FZ(I)*Y0Z
+ TPX= TX(I)*X0X + TY(I)*X0Y + TZ(I)*X0Z
+ TPY= TX(I)*Y0X + TY(I)*Y0Y + TZ(I)*Y0Z
+C To avoid rounding error in looking for expansion limits
+ DUM = ABS(FPY)
+ DUM1= ABS(TPY)
+ IF(DUM.LT.R4PREC) FPY=0.0
+ IF(DUM1.LT.R4PREC) TPY=0.0
+ YMAX= DUM+DUM1
+ YMIN=-YMAX
+C ################ non serve loppando sulle celle e non sy yesp?##
+C To avoid yesp value exactly on grid boundaries, which can cause
+C the skipping of some cells. If on cell border is mouved a bit
+C into the cell. (If expanded over the previous cell it doesn't
+C find y intersection with projected surf.elem skipping this y cell)
+C IF(MOD(YMIN+YY0,DCELL).LT.(DCELL/1000.)) YMIN=YMIN+(DCELL/1000.)
+C YMAX can't be expanded to the cell boundary . In the last cell
+C YESP must be <=YMAX or no y intersection will be found.
+C
+C First and last Y cell covered by the surf. element, in which
+C YMAX and YMIN falls; we have a little asimmetry here:
+C if Ymin is on left cell boundary this doesn't expand
+C Y on the previous cell. But if Ymax is on the rigth cell
+C boundary then expand Y over the next cell.
+ JYMAX=INT((YMAX+YY0)/DCELL)+1
+ JYMIN=INT((YMIN+YY0)/DCELL)+1
+C Test if these two cells are on the projection plane boundaries
+ IF(JYMIN.LE.0) JYMIN=1
+ IF(JYMAX.GT.NCELL) THEN
+ JYMAX=NCELL
+ YMAX=CELLMAX-CELLMIN-YY0
+ ENDIF
+C
+C ------------------- loop on y expansion grid range
+C
+ DO 80 JY=JYMIN,JYMAX
+C An Y value YESP is computed for each cell covered by the surf.el.
+C then the intersections between the projected surf.el. boundary
+C and the line Y=YESP are found: these are the X expansion limits
+C for this JY row of proj. plane.
+C Potrebbe essere utile mettere YESP vicino al bordo della
+C cella piu' vicino al punto YY0 in cui cade l'elemento di
+C superficie, cioe' al bordo superiore sotto YY0 e a quello
+C inferiore sopra. In questo modo ci si mette nelle condizioni
+C di trovare le intersezioni X piu' distanti possibile fra loro
+C e quindi di espandersi sul massimo numero di celle, qui invece
+C puo' capitare di avere YESP che supera l'espasione del surf. el.
+C e quindi non trova intersezioni anche se in effetti il surf. el.
+C entra un poco nella cella.
+C
+ IF(JY.NE.JYMAX) THEN
+ YESP=YMIN+(JY-JYMIN)*DCELL
+ ELSE
+C If YESP is out of the plane boundary it's set on boundary
+ YESP=YMAX
+ ENDIF
+C
+C Intesecting the line Y=yesp with the boundary of the proj.surf. el.
+ IXESP=1
+ IF(TPY.NE.0.) THEN
+ H=(YESP-FPY)/TPY
+ IF(ABS(H).LE.1.) THEN
+ XESP12(IXESP)=FPX+TPX*H
+ IXESP=IXESP+1
+ ENDIF
+ H=(YESP+FPY)/TPY
+ IF(ABS(H).LE.1.) THEN
+ XESP12(IXESP)=-FPX+TPX*H
+ IXESP=IXESP+1
+ ENDIF
+ ENDIF
+ IF(FPY.NE.0.) THEN
+ H=(YESP-TPY)/FPY
+ IF(ABS(H).LE.1.) THEN
+ XESP12(IXESP)=TPX+FPX*H
+ IXESP=IXESP+1
+ ENDIF
+ H=(YESP+TPY)/FPY
+ IF(ABS(H).LE.1.) THEN
+ XESP12(IXESP)=-TPX+FPX*H
+ IXESP=IXESP+1
+ ENDIF
+ ENDIF
+C
+ IF(IXESP.EQ.1) THEN
+C No intersection, this can happen if the surf element is
+C all into Dcell, this case is considered at the end of
+C the expansion loop. This can also happen at the Y proj.surf.el.
+C boundary if YESP (about on cell boundary) falls is just outside
+C the cell it should represent.
+ GOTO 80
+ ELSE IF(IXESP.EQ.2) THEN
+C Only one intersection, typically at y surf.el. boundary, when both
+C are at the edge of the surf.elem. In this case we should have
+C two intersection of equal value, but rounding errors can alter.
+ XMAX=XESP1
+ XMIN=XESP1
+ XESP2=XESP1
+ ELSE IF(IXESP.GT.3) THEN
+C Four intersection are too mutch; this can happen if YESP falls
+C EXACTELY on the proj.surf.el. diagonal, or when TPY or FPY are
+C about 0, in this case rounding errors can give spurious intersections
+C with H and YESP about 0.
+ WRITE(6,1997) IF,I,IXESP-1,(XESP12(J),J=1,IXESP-1)
+ IF(PRINTFILE) WRITE(NFILE,1997)
+ 1 IF,I,IXESP-1,(XESP12(J),J=1,IXESP-1)
+ 1997 FORMAT(' WARINIG! : in light calculation, phase:',I4,' point:'
+ 1 I6/' Too mutch expansion directions:',I3,1X,4E15.5)
+ IF(IXESP.EQ.4) THEN
+ XMAX=MAX(XESP1,XESP2,XESP3)
+ XMIN=MAX(XESP1,XESP2,XESP3)
+ ELSE
+ XMAX=MAX(XESP1,XESP2,XESP3,XESP4)
+ XMIN=MAX(XESP1,XESP2,XESP3,XESP4)
+ ENDIF
+ ELSE
+C IXESP=3 : two intersections, wich are the X expansion limits
+ XMAX=MAX(XESP1,XESP2)
+ XMIN=MIN(XESP1,XESP2)
+ ENDIF
+C
+ JXMIN=INT((XMIN+XX0)/DCELL)+1
+ JXMAX=INT((XMAX+XX0)/DCELL)+1
+ IF(JXMIN.LE.0) JXMIN=1
+ IF(JXMAX.GT.NCELL) JXMAX=NCELL
+C
+C ------------------- loop on X expansion grid range
+ DO 82 JX=JXMIN,JXMAX
+C
+ J=ICELL(JX,JY)
+ NCOUNTE=NCOUNTE+1
+C
+C If the grid point JX,JY is void then expand the I point over it
+ IF(J.EQ.0) THEN
+ ICELL(JX,JY)=I
+ IGAINED(I)=IGAINED(I)+1
+ NCOUNTV=NCOUNTV+1
+C If filled by a different object
+ ELSE IF(NUMOBJ(J).NE.NUMOBJ(I)) THEN
+ NCOUNTD=NCOUNTD+1
+ DISTI=X(I)*Z0X+Y(I)*Z0Y+Z(I)*Z0Z
+ DISTJ=X(J)*Z0X+Y(J)*Z0Y+Z(J)*Z0Z
+C if this object is the nearer (the nearer has the greatest DIST)
+ IF(DISTI.GT.DISTJ) THEN
+ ICELL(JX,JY)=I
+ IGAINED(I)=IGAINED(I)+1
+ ILOST(J)=ILOST(J)+1
+ NCOUNTD1=NCOUNTD1+1
+ ELSE
+C The other is the nearer
+ IGAINED(I)=IGAINED(I)+1
+ ILOST(I)=ILOST(I)+1
+ NCOUNTD2=NCOUNTD2+1
+ ENDIF
+ ELSE
+C filled by the same object
+C the cell remains assigned to the other object
+ NCOUNTF=NCOUNTF+1
+ ENDIF
+C
+ 82 CONTINUE
+C ------------------------------- END OF X EXPANSION LOOP 82
+ 80 CONTINUE
+C ------------------------------- END OF Y EXPANSION LOOP 80
+C
+C Test if this point has some cell which is not eclipsed and in which
+C there aren't other points of the same object
+ IF (IGAINED(I).GT.ILOST(I)) GOTO 30
+C K is the surface element of the cell in which this point falls
+ K=ICELL(IX,IY)
+ IF(K.LE.0) THEN
+C This is the case in which the surf.elem. << dcell has been jumped
+C by the loop. This can hapen when no interseciotn is found, being
+C the surf.el. << dcell.
+ ICELL(IX,IY)=I
+ IGAINED(I)=IGAINED(I)+1
+ NCOUNTV=NCOUNTV+1
+C If filled by a different object (could be if surf.el.<<dcell, as before)
+C But generally this has been done before, we have redundancy here.
+ ELSE IF(NUMOBJ(K).NE.NUMOBJ(I)) THEN
+ NCOUNTD=NCOUNTD+1
+ DISTI=X(I)*Z0X+Y(I)*Z0Y+Z(I)*Z0Z
+ DISTK=X(K)*Z0X+Y(K)*Z0Y+Z(K)*Z0Z
+C if this object is the nearer (the nearer has the greatest DIST)
+ IF(DISTI.GT.DISTK) THEN
+ ICELL(IX,IY)=I
+ IGAINED(I)=IGAINED(I)+1
+ ILOST(K)=ILOST(K)+1
+ NCOUNTD1=NCOUNTD1+1
+ ELSE
+C The other is the nearer
+ IGAINED(I)=IGAINED(I)+1
+ ILOST(I)=ILOST(I)+1
+ NCOUNTD2=NCOUNTD2+1
+ ENDIF
+ ELSE
+C Else: last chance for this point: its cell is assigned to
+C a surf. el. of the same object, it enters into its linklist
+ NEXT(I)=NEXT(K)
+ NEXT(K)=I
+ NCOUNTL=NCOUNTL+1
+ ENDIF
+C
+ 30 CONTINUE
+C ---------------------------------- End of surf. elem. loop
+ IF(LOST.GT.0) THEN
+ WRITE(6,1000) LOST
+ IF(PRINTFILE) WRITE(NFILE,1000) LOST
+ 1000 FORMAT(' ERROR: ',I6,' LIGHT POINTS LOST IN PROJECTION '/
+ 1 ' This is a real program bug, because you should never have'/
+ 2 ' caused this error in subroutine LUCE')
+ ENDIF
+C .........................Finally computes the ligth to the
+C observer by integration over the grid
+C
+C ---------------------------------- BEGIN LIGHT SUM LOOP 40
+ DO 40 I=1,NPNT
+ IF(IGAINED(I).LE.0) GOTO 40
+ IF(NEXT(I).GT.0) NCOUNTL1=NCOUNTL1+1
+C not eclipsed fraction of surf.element
+ REMAINED=1.- REAL(ILOST(I))/REAL(IGAINED(I))
+ IF(REMAINED.LE.0.0) GOTO 40
+C ----------- loop on the link list of point I
+C The eclipsed fraction of point I is used for the whole link list
+ K=I
+ 45 CONTINUE
+ GEOMETRY=(1.0+XLIMB(K)*(COSG(K)-1.0))*COSG(K) * REMAINED
+ ALUCE= GEOMETRY *(FB(K)+FRIFL(K))
+ ALUCER=GEOMETRY * FRIFL(K)
+C A concave object is represented by a number of convex parts,
+C each drawn as a different object and numbered: obj.num ,obj.num+3 ,..+6...
+C KOBJ= 1 for first,2 for second, 3 for third object.
+ KOBJ=NUMOBJ(K)
+C The following holds also for numobj<=3, but if instead of -/* saves time
+ IF(KOBJ.GT.3) KOBJ=KOBJ-((KOBJ-1)/3)*3
+ CBOL(KOBJ,IF) =CBOL(KOBJ,IF) +ALUCE
+ CBOLR(KOBJ,IF)=CBOLR(KOBJ,IF)+ALUCER
+ DO 50 IB=1,NBANDE
+ IF(IBAND(IB).LE.0) GOTO 50
+ ALUCE=GEOMETRY * F(K,IB)
+ C(KOBJ,IB,IF)=C(KOBJ,IB,IF)+ALUCE
+ 50 CONTINUE
+ NCOUNTT=NCOUNTT+1
+ IF(NEXT(K).GT.0) THEN
+ K=NEXT(K)
+ NCOUNTTL=NCOUNTTL+1
+ GOTO 45
+ ENDIF
+C --------------- end of loop on the link list
+ 40 CONTINUE
+C ------------------------------ END OF THE LIGHT LOOP
+C
+ CBOLT(IF)=CBOL(1,IF)+CBOL(2,IF)+CBOL(3,IF)
+ CBOLTR(IF)=CBOLR(1,IF)+CBOLR(2,IF)+CBOLR(3,IF)
+ DO 60 IB=1,NBANDE
+ IF(IBAND(IB).GT.0)
+ 1 CT(IB,IF)=C(1,IB,IF)+C(2,IB,IF)+C(3,IB,IF)
+ 60 CONTINUE
+C
+C ---------------------------------------------- PRINTING GRID
+C
+C ............ printing the grid points at each nprint phases
+C .................. Printing the object in each grid
+ IF(NPRINT.LE.0.OR..NOT.PRINTFILE) GOTO 100
+ IF(MOD(IF,NPRINT).EQ.0) THEN
+C ............ printing statistical data
+ IF(NPRINTL.LT.1) GOTO 100
+ WRITE(NFILE,5900) IF,NCOUNTS,NCOUNTL,NCOUNTS-NCOUNTL,
+ 1 NCOUNTE,NCOUNTV,NCOUNTD,NCOUNTD1,NCOUNTD2,
+ 2 NCOUNTF,
+ 3 NCOUNTT,NCOUNTTL,NCOUNTT-NCOUNTTL,NCOUNTL1
+ 5900 FORMAT(/' Statistics summary for ligth computation in phase:',I5/
+ 1/' Number of considered surface elements:',I10/
+ 2' of which: after the first in link lists:',I10/
+ 2' of which: not in link lists or first of a list:',I10/
+ 4' Number of cell expansions:',I10/
+ 4' Number of expansions over void cells:',I10/
+ 4' Number of expansions over different object''s cells:',I10/
+ 4' of which resulting in other object''s eclipsing :',I10/
+ 4' of which resulting in this object''s eclipsing :',I10/
+ 4' Number of expansions on same object''s cells:',I10/
+ 9' Total number of remaining surface elements:',I10/
+ 1' of which:in linklists :',I10/
+ 1' of which:not in linklists :',I10/
+ 1' Total number of remaining linklists :',I10)
+C
+ IF(NPRINTL.LT.2) GOTO 100
+ WRITE(NFILE,5999) IF,FASE(IF),FASE(IF)/6.28318
+ 5999 FORMAT(/' Number of object in each grid point for fase n.',I5,
+ 1 ' => ',2G15.7)
+ NPSTP=NCELL/100
+ NPINIT=1
+ DO 62 IXX=1,NPSTP
+ NPFIN=NPINIT+(100-1)
+ DO 63 IY=1,NCELL
+ WRITE(NFILE,6000)IY,YCELL(IY),
+ 1 (NUMOBJ(ABS(ICELL(IX,IY))),IX=NPINIT,NPFIN)
+ 6000 FORMAT(1X,I3,1X,G12.3,' = ',100I1)
+ 63 CONTINUE
+ NPINIT=NPFIN+1
+ 62 CONTINUE
+ IF(NPINIT.LE.NCELL) THEN
+ DO 64 IY=1,NCELL
+ WRITE(NFILE,6000)IY,YCELL(IY),
+ 1 (NUMOBJ(ABS(ICELL(IX,IY))),IX=NPINIT,NCELL)
+ 64 CONTINUE
+ ENDIF
+C .................. Prints surf. point for each grid point
+ IF(NPRINTL.LT.3) GOTO 100
+C Format I5 can't contain numbers > 99999 ,in this case it doesn't print
+ IF(NPNT.GT.99999) GOTO 100
+ WRITE(NFILE,6001) IF,FASE(IF),FASE(IF)/6.283185
+ 6001 FORMAT(' Number of surf.element in each grid point for fase n.',
+ 1 I5,' => ',2G15.7)
+ NPSTP=NCELL/20
+ NPINIT=1
+ DO 67 IXX=1,NPSTP
+ NPFIN=NPINIT+(20-1)
+ WRITE(NFILE,6004)
+ WRITE(NFILE,6002)(IX,IX=NPINIT,NPFIN)
+ 6002 FORMAT(' Cell grid numbers: ',20(1X,I4))
+ WRITE(NFILE,6003)(XCELL(IX),IX=NPINIT,NPFIN)
+ 6003 FORMAT(' Cell grid center :',20(F5.2))
+ WRITE(NFILE,6004)
+ 6004 FORMAT(1X,119(1H-))
+ DO 68 IY=1,NCELL
+ WRITE(NFILE,6005) IY,YCELL(IY),
+ 1 (ICELL(IX,IY),IX=NPINIT,NPFIN)
+ 6005 FORMAT(1X,I3,1X,G12.3,' = ',20I5)
+ 68 CONTINUE
+ NPINIT=NPFIN+1
+ 67 CONTINUE
+ IF(NPINIT.LE.NCELL) THEN
+ WRITE(NFILE,6004)
+ WRITE(NFILE,6002)(IX,IX=NPINIT,NCELL)
+ WRITE(NFILE,6003)(XCELL(IX),IX=NPINIT,NCELL)
+ WRITE(NFILE,6004)
+ DO 69 IY=1,NCELL
+ WRITE(NFILE,6005) IY,YCELL(IY),
+ 1 (ICELL(IX,IY),IX=NPINIT,NCELL)
+ 69 CONTINUE
+ ENDIF
+ WRITE(NFILE,6004)
+C
+C ........................... Printing eclipse-suppressed points
+C ############################# solo utili per test non devono esistere
+ IF(NPRINTL.LT.4) GOTO 100
+ WRITE(NFILE,6010)
+ 6010 FORMAT(' Surface points which appear in the grid, but aren''t',
+ 1 ' considered, being eclipsed:')
+ K=0
+ NDEL=0
+ DO 90 IY=1,NCELL
+ DO 90 IX=1,NCELL
+ KK=ICELL(IX,IY)
+ IF(KK.LE.0) GOTO 90
+ IF(IGAINED(KK).GT.ILOST(KK)) GOTO 90
+ NDEL=NDEL+1
+ K=K+1
+ NBUF(K)=KK
+ IF(K.EQ.MAXBUF) THEN
+ WRITE(NFILE,6015) (NBUF(J),J=1,MAXBUF)
+ 6015 FORMAT(1X,20I5)
+ K=0
+ ENDIF
+ 90 CONTINUE
+ IF(K.GT.0) WRITE(NFILE,6015) (NBUF(J),J=1,K)
+ WRITE(NFILE,6017) NDEL
+ 6017 FORMAT(' total:',I5,' surf.points not considered in the grid')
+C
+C ............................... Printing the link lists
+C
+ IF(NPRINTL.LT.5) GOTO 100
+ WRITE(NFILE,6020)
+ 6020 FORMAT(' Linklists for the surf.points in the grid:')
+ K=0
+ NLIST=0
+ DO 95 I=1,NPNT
+ KK=I
+ IF(KK.LE.0) GOTO 95
+ IF(IGAINED(KK).LE.0) GOTO 95
+ IF(IGAINED(KK).LE.ILOST(KK) ) GOTO 95
+ IF(NEXT(KK).LE.0) GOTO 95
+C A linklist has been found, loops on its points.
+ NLIST=NLIST+1
+ WRITE(NFILE,6025) NLIST,KK
+ 6025 FORMAT(' Linklist:',I5,' First surf.element:',I6)
+ 96 K=K+1
+ NBUF(K)=KK
+ IF(K.EQ.MAXBUF) THEN
+ WRITE(NFILE,6015) (NBUF(J),J=1,MAXBUF)
+ K=0
+ ENDIF
+ KK=NEXT(KK)
+ IF(KK.GT.0) GOTO 96
+ IF(K.GT.0) THEN
+ WRITE(NFILE,6015) (NBUF(J),J=1,K)
+ K=0
+ ENDIF
+ 95 CONTINUE
+ IF(NLIST.LE.0) WRITE(NFILE,6030)
+ 6030 FORMAT(' NONE')
+C
+ ENDIF
+ 100 CONTINUE
+C ---------------------------------------------- END OF PRINTING
+C
+ 10 CONTINUE
+C ................................. end of loop on phases
+C
+ IF(LOST.GT.0) THEN
+ WRITE(6,1100) LOST
+ IF(PRINTFILE) WRITE(NFILE,1100) LOST
+ 1100 FORMAT(' ERROR: ',I6,' LIGHT POINTS LOST IN INTEGRATION '/
+ 1 ' This is a real program bug, because you should never have'/
+ 2 ' caused this error in subroutine LUCE')
+ ENDIF
+C
+ RETURN
+ END
+C
+ SUBROUTINE MAXABS(XMAX,N,X)
+C --------------------------------------------------------
+C returns xmax, maximum absolute value of x(n)
+C --------------------------------------------------------
+ DIMENSION X(N)
+ XMAX=ABS(X(1))
+ DO 10 I=2,N
+ XX=X(I)
+ IF(XX.LT.0.) XX=-XX
+ IF(XX.GT.XMAX) XMAX=XX
+ 10 CONTINUE
+ RETURN
+ END
+C
+C
+ SUBROUTINE MAXMIN(XMAX,XMIN,N,X)
+C -----------------------------------------------------------
+C This routine find the maximum and minimum values for X
+C ROUTINE NOT USED!
+C -----------------------------------------------------------
+ DIMENSION X(N)
+ XMAX=X(1)
+ XMIN=X(1)
+ DO 10 I=2,N
+ IF(XMAX.LT.X(I)) XMAX=X(I)
+ IF(XMIN.GT.X(I)) XMIN=X(I)
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE NOREFL(NPNT,NIF,FBOLREF,FRIFL)
+C --------------------------------------------------
+C This routine sets the reflected flux to zero
+C REALLY, MAY BE SIMPLIFIED !!!
+C ------------------------------------------------
+ DIMENSION FRIFL(NPNT),FBOLREF(3),NIF(2,3)
+C
+C ..................... loop on objects
+ DO 10 I=1,3
+ N1=NIF(1,I)
+ N2=NIF(2,I)
+ FBOLREF(I)=0.0
+ IF(N1.LE.0.OR.N2.LE.0) GOTO 10
+ DO 20 J=N1,N2
+ FRIFL(J)=0.0
+ 20 CONTINUE
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE NORM(N,A,B)
+C --------------------------
+C Multiplies B(i) by A
+C --------------------------
+ DIMENSION B(N)
+ DO 10 I=1,N
+ 10 B(I)=B(I)*A
+ RETURN
+ END
+C
+ SUBROUTINE NORMC(ANORMC,IFLAG)
+C ---------------------------------------------------------------
+C Normalizes the light curve to have max value= ANORMC
+C This is done for bolometric flux CBOLT,
+C and each requested color band CT,(requested bands have iflag>0)
+C ---------------------------------------------------------------
+ PARAMETER (MAXFASI=1000, MAXBANDE=5)
+ COMMON /LIGHT/ NFASIMX,NFASI,NBANDE,
+ 1 FASE(MAXFASI),O(MAXBANDE,MAXFASI),
+ 2 C(3,MAXBANDE,MAXFASI),CT(MAXBANDE,MAXFASI),
+ 3 CBOL(3,MAXFASI),CBOLT(MAXFASI),
+ 4 CBOLR(3,MAXFASI),CBOLTR(MAXFASI),
+ 5 AREA(3),FBOLTOT(3),FBOL(3),FBOLREF(3)
+C
+ DIMENSION IFLAG(MAXBANDE)
+C
+C normalizing bolometric flux, total and for each object
+C also reflection sources are normalized
+ AMAX=AMAX2(NFASI,CBOLT)
+ IF(AMAX.EQ.0.0) GOTO 100
+ DO 10 I=1,NFASI
+ CBOLT(I)=CBOLT(I)/AMAX*ANORMC
+ CBOL(1,I)=CBOL(1,I)/AMAX*ANORMC
+ CBOL(2,I)=CBOL(2,I)/AMAX*ANORMC
+ CBOL(3,I)=CBOL(3,I)/AMAX*ANORMC
+ CBOLTR(I)=CBOLTR(I)/AMAX*ANORMC
+ CBOLR(1,I)=CBOLR(1,I)/AMAX*ANORMC
+ CBOLR(2,I)=CBOLR(2,I)/AMAX*ANORMC
+ CBOLR(3,I)=CBOLR(3,I)/AMAX*ANORMC
+ 10 CONTINUE
+ 100 CONTINUE
+C
+C normalizing color band flux, total and for each object
+ DO 20 I=1,NBANDE
+ IF(IFLAG(I).LE.0) GOTO 20
+ AMAXB=CT(I,1)
+ DO 30 J=2,NFASI
+ IF(CT(I,J).GT.AMAXB) AMAXB=CT(I,J)
+ 30 CONTINUE
+ IF(AMAXB.EQ.0.0) GOTO 20
+ DO 40 J=1,NFASI
+ CT(I,J)=CT(I,J)/AMAXB*ANORMC
+ C(1,I,J)=C(1,I,J)/AMAXB*ANORMC
+ C(2,I,J)=C(2,I,J)/AMAXB*ANORMC
+ C(3,I,J)=C(3,I,J)/AMAXB*ANORMC
+ 40 CONTINUE
+ 20 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE NORML(NPNT,FINPUT,TL,F)
+C -----------------------------------------------------------
+C Normalization routine. Normalizes the total flux to an
+C input given value FINPUT. ( Here the flux of each surface
+C element had been previously multiplied by the element area)
+C -----------------------------------------------------------
+ DIMENSION F(NPNT)
+ TL=0.0
+ DO 10 I=1,NPNT
+ 10 TL=TL+F(I)
+ DO 20 I=1,NPNT
+ 20 F(I)=(FINPUT/TL)*F(I)
+ RETURN
+ END
+C
+ SUBROUTINE NULL(N,A)
+C -----------------------
+C Fills A with 0.0
+C ----------------------
+ DIMENSION A(N)
+ DO 10 I=1,N
+ 10 A(I)=0.0
+ RETURN
+ END
+C
+ FUNCTION PLANK(AL,TEMPER)
+C -------------------------------------------------
+C Value of Plank function, for lambda AL ( Anstrong )
+C and temperature TEMPER ( K )
+C PLANK= erg/cm**2/sec/sterad / lambda unit
+C -------------------------------------------------
+C
+C This is for lambda in Anstrong : 1 cm = 1.E8 Anstrong, may give overflow
+C PLANK=1.1909D35 /
+C 1 ( AL**5 * (DEXP (1.4398D+8/(DBLE(AL)*DBLE(TEMPER)) ) -1.D0) )
+C This is for lambda in cm
+C PLANK=1.1909D-5 /
+C 1 ( AL**5 * (DEXP (1.4398D0/(DBLE(AL)*DBLE(TEMPER)) ) -1.D0) )
+C This is for lambda in A/1000 and A
+ AL1=AL/1000.
+ PLANK=1.1909D20 /
+ 1 ( AL1**5 * (DEXP (1.4398D+8/(DBLE(AL)*DBLE(TEMPER)) ) -1.D0) )
+ RETURN
+ END
+C
+ FUNCTION PLANKL(AL,TEMPER)
+C -------------------------------------------------
+C Value of Plank function * delta lambda
+C delta lambda =1 A
+C -------------------------------------------------
+C
+C This is for lambda in Anstrong : 1 cm = 1.E8 Anstrong
+ AL1=AL/1000.
+ PLANKL=1.1909D12 /
+ 1 ( AL1**5 * (DEXP (1.4398D+8/(DBLE(AL)*DBLE(TEMPER)) ) -1.D0) )
+C PLANKL=1.1909D27 /
+C 1 ( AL**5 * (DEXP (1.4398D+8/(DBLE(AL)*DBLE(TEMPER)) ) -1.D0) )
+ RETURN
+ END
+C
+ SUBROUTINE PLOTTING(PRINTFILE,NFILE,PRINT6,
+ 1 AMOUNTP,NFILE6P,NFILEP,NASCIIP)
+C ----------------------------------------------------------------
+C Plotting routine : now only writes on unit nasciip,
+C plots are done by an'other program.
+C ----------------------------------------------------------------
+ PARAMETER (MAXFASI=1000, MAXBANDE=5)
+ COMMON /LIGHT/ NFASIMX,NFASI,NBANDE,
+ 1 FASE(MAXFASI),O(MAXBANDE,MAXFASI),
+ 2 C(3,MAXBANDE,MAXFASI),CT(MAXBANDE,MAXFASI),
+ 3 CBOL(3,MAXFASI),CBOLT(MAXFASI),
+ 4 CBOLR(3,MAXFASI),CBOLTR(MAXFASI),
+ 5 AREA(3),FBOLTOT(3),FBOL(3),FBOLREF(3)
+C
+ LOGICAL PRINTFILE,PRINT6
+ DIMENSION FASI01(MAXFASI)
+C
+ DO 10 I=1,NFASI
+ 10 FASI01(I)=FASE(I)/(2.*3.141593)
+C
+ IF(NASCIIP.GT.0.AND.NASCIIP.LE.99) THEN
+ OPEN(UNIT=NASCIIP,FORM='FORMATTED',STATUS='NEW')
+ DO 20 I=1,NFASI
+ WRITE(NASCIIP,1000) I,FASI01(I),CBOLT(I),
+ 1 (CT(J,I),J=1,NBANDE),(O(J1,I),J1=1,NBANDE)
+ 1000 FORMAT(1X,I5,1X,F6.4,11(1X,F9.6))
+ 20 CONTINUE
+ CLOSE(UNIT=NASCIIP)
+ IF(PRINT6) WRITE(6,2000) NASCIIP,NFASI,NBANDE
+ IF(PRINTFILE) WRITE(NFILE,2000) NASCIIP,NFASI,NBANDE
+ 2000 FORMAT(/' Light curve written on unit:',I2,' Num.phases:',
+ 1 I8,' Num. of colors:',I2)
+ ENDIF
+C
+ RETURN
+ END
+C
+ SUBROUTINE POINT(PRINT6,PRINTFILE,NFILE,NRF,NTHF,PKRI,NOBJ,
+ 1 NPUNT,NPUNTC,R,RY,X0,Y0,Z0,ALZO,ROTAZ,R4PREC,
+ 2 NPNTMX,X,Y,Z,G,GX,GY,GZ,FX,FY,FZ,TX,TY,TZ,AREA,A,NUMCOARSE,
+ 3 NUMOBJ)
+C ----------------------------------------------------------------
+C Draws a bright point surface, in a rectangle
+C R=x half-edge
+C RY=y half-edge
+C NRF= number of fine r mesh in which 2*R is divided.
+C NTHF= coarsing factor=number of fine for each coarse
+C NRC= number of r coarse
+C if PKRI > 0.0 also the back surface of the rectangle is drawn.
+C
+C After drawing on the plane x-y the cylinder can be rotated:
+C ROTAZ= anti-clockwise angle of rotazion on the plane (x-x'angle)
+C ALZO= clockwise rotation around the y' axis (z-z" angle)
+C
+C ----------------------------------------------------------------
+ DIMENSION X(NPNTMX),Y(NPNTMX),Z(NPNTMX)
+ DIMENSION G(NPNTMX),GX(NPNTMX),GY(NPNTMX),GZ(NPNTMX)
+ DIMENSION FX(NPNTMX),FY(NPNTMX),FZ(NPNTMX)
+ DIMENSION TX(NPNTMX),TY(NPNTMX),TZ(NPNTMX)
+ DIMENSION A(NPNTMX)
+ DIMENSION NUMCOARSE(NPNTMX),NUMOBJ(NPNTMX)
+ LOGICAL PRINTFILE,PRINT6
+C
+C Coarse mesh
+ NRC=NRF/NTHF
+ IF(MOD(NRF,NTHF).NE.0) THEN
+ WRITE(6,1000) NRF,NTHF
+ 1000 FORMAT(' ERROR! In drawing a rectangle the number of r mesh:',
+ 1 I5,' is not consistent with the coarsing factor:',I5)
+ ENDIF
+C
+C coarse mesh
+ EDGE=2.0*R
+ DRC=EDGE/NRC
+ DRC5=DRC*0.5
+C fine r mesh
+ DRF=DRC/NTHF
+ DRF5=DRF*0.5
+C
+C NPUNTC will be incremented below. (NPUNTC is the first void coarse)
+C IC is instead the actual coarse, or the last filled one.
+C NPUNT is the number of surf.el. generated by this routine
+ IF(PKRI.GT.0.0) THEN
+ IC=NPUNTC-2
+ ELSE
+ IC=NPUNTC-1
+ ENDIF
+C
+C Y coarse and fine mesh (Y must be a multiple of X: DRCY=DRC)
+ DRCY=DRC
+ EDGEY=2.0*RY
+ NRCY=EDGEY/DRCY
+ DRCY5=DRCY*0.5
+ DRFY=DRCY/NTHF
+ DRFY5=DRFY*0.5
+C ...............
+ NPUNT=0
+ AREA=0.0
+C
+C ................... Coarse X loop ( from center -R to R)
+ DO 10 IXC=1,NRC
+ X1C=-R+(IXC-1)*DRC
+ X2C= X1C+DRC
+C
+C ................... Coarse Y loop ( from center -R to R)
+ DO 10 IYC=1,NRCY
+ Y1C=-RY+(IYC-1)*DRCY
+ Y2C= Y1C+DRCY
+C
+ IF(PKRI.GT.0.0) THEN
+ IC=IC+2
+ ELSE
+ IC=IC+1
+ ENDIF
+C
+C ................... Fine X loop ( from X1C to X2C )
+ DO 10 IXF=1,NTHF
+ XF=X1C+(IXF-1)* DRF +DRF5
+C
+C ................... Fine Y loop ( from Y1C to Y2C )
+ DO 10 IYF=1,NTHF
+ YF=Y1C+(IYF-1)* DRFY +DRFY5
+C
+C Upper surface point
+ IF(NPUNT.GE.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 100
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)= XF
+ Y(NPUNT)= YF
+ Z(NPUNT)= 0.0
+ GX(NPUNT)= 0.0
+ GY(NPUNT)= 0.0
+ GZ(NPUNT)= 1.0
+ FX(NPUNT)= 0.
+ FY(NPUNT)= DRFY5
+ FZ(NPUNT)= 0.
+ TX(NPUNT)= DRF5
+ TY(NPUNT)= 0.
+ TZ(NPUNT)= 0.
+ A(NPUNT) = DRF*DRFY
+ NUMCOARSE(NPUNT)=IC
+ NUMOBJ(NPUNT)=NOBJ
+C area computation
+ AREA=A(NPUNT)+AREA
+C
+C Lower surface point
+ 100 IF(PKRI.LE.0.0) GOTO 10
+ IF(NPUNT.GE.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 10
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)= XF
+ Y(NPUNT)= YF
+ Z(NPUNT)= 0.0
+ GX(NPUNT)= 0.0
+ GY(NPUNT)= 0.0
+ GZ(NPUNT)= -1.0
+ FX(NPUNT)= 0.0
+ FY(NPUNT)= DRFY5
+ FZ(NPUNT)= 0.0
+ TX(NPUNT)= -DRF5
+ TY(NPUNT)= 0.0
+ TZ(NPUNT)= 0.0
+ A(NPUNT) = DRF*DRFY
+ NUMCOARSE(NPUNT)=IC+1
+ NUMOBJ(NPUNT)=NOBJ
+C Total area computation
+ AREA=A(NPUNT)+AREA
+C
+ 10 CONTINUE
+C
+C IF the IC+1 coarse has been filled by loop 10
+ IF(PKRI.GT.0.0) IC=IC+1
+C
+ NPUNTC=IC+1
+C
+C Constant potential USEFUL FOR POINT ?
+ DO 59 I=1,NPUNT
+ 59 G(I)=1.
+C
+ IF(LOST.GT.0) THEN
+ WRITE(6,2000) NPUNT,LOST
+ IF(PRINTFILE) WRITE(NFILE,1000) NPUNT,LOST
+ 2000 FORMAT(' ERROR! In POINT:OBJECT NOT PROPERLY DRAWN! '/
+ 1 ' space available for only:',I6,' surface points.',
+ 2 ' lost points:',I6/
+ 3 ' REDUCE THE MESH PARAMETER OR GO INTO THE FORTRAN LIST TO',
+ 4 ' INCREASE MAXPT')
+ ENDIF
+C
+C ................ optional rotation .....
+C
+ IF(ALZO.EQ.0.0.AND.ROTAZ.EQ.0.0) GOTO 800
+C The disk is rotated by rotaz around the z axis (anti-clockwise)
+C then around the new Y axis (clockwise) by alzo
+ COSP=COS(ROTAZ)
+ SINP=SIN(ROTAZ)
+ COSI=COS(ALZO)
+ SINI=SIN(ALZO)
+C
+C To prevent precision broblems
+ IF(ABS(COSP).LT.R4PREC) COSP=0.0
+ IF(ABS(SINP).LT.R4PREC) SINP=0.0
+ IF(ABS(COSI).LT.R4PREC) COSI=0.0
+ IF(ABS(SINI).LT.R4PREC) SINI=0.0
+C
+ DO 80 I=1,NPUNT
+ X1= X(I)* (COSP*COSI) + Y(I) * SINP - Z(I) * (COSP*SINI)
+ Y1=-X(I)* (SINP*COSI) + Y(I) * COSP + Z(I) * (SINP*SINI)
+ Z1= X(I)* SINI + Z(I) * COSI
+ X(I)=X1
+ Y(I)=Y1
+ Z(I)=Z1
+ X1= GX(I)* (COSP*COSI) + GY(I) * SINP - GZ(I) * (COSP*SINI)
+ Y1=-GX(I)* (SINP*COSI) + GY(I) * COSP + GZ(I) * (SINP*SINI)
+ Z1= GX(I)* SINI + GZ(I) * COSI
+ GX(I)=X1
+ GY(I)=Y1
+ GZ(I)=Z1
+ X1= TX(I)* (COSP*COSI) + TY(I) * SINP - TZ(I) * (COSP*SINI)
+ Y1=-TX(I)* (SINP*COSI) + TY(I) * COSP + TZ(I) * (SINP*SINI)
+ Z1= TX(I)* SINI + TZ(I) * COSI
+ TX(I)=X1
+ TY(I)=Y1
+ TZ(I)=Z1
+ X1= FX(I)* (COSP*COSI) + FY(I) * SINP - FZ(I) * (COSP*SINI)
+ Y1=-FX(I)* (SINP*COSI) + FY(I) * COSP + FZ(I) * (SINP*SINI)
+ Z1= FX(I)* SINI + FZ(I) * COSI
+ FX(I)=X1
+ FY(I)=Y1
+ FZ(I)=Z1
+C
+ 80 CONTINUE
+C
+ 800 CONTINUE
+C shifting the point from the origin=(0,0,0) to x0,y0,z0
+ IF(X0.NE.0.0) THEN
+ DO 90 I=1,NPUNT
+ 90 X(I)=X(I)+X0
+ ENDIF
+ IF(Y0.NE.0.0) THEN
+ DO 91 I=1,NPUNT
+ 91 Y(I)=Y(I)+Y0
+ ENDIF
+ IF(Z0.NE.0.0) THEN
+ DO 92 I=1,NPUNT
+ 92 Z(I)=Z(I)+Z0
+ ENDIF
+C
+ RETURN
+ END
+C
+ SUBROUTINE PRINTING(AMOUNT,NT6,NT,IPPAR1,IPPAR2,
+ 1 NAMFLAG,DESCFLAG,NAMPAR,DESCPAR,AUS)
+C ----------------------------------
+C Output printing routine
+C ----------------------------------
+ PARAMETER (MAXPT=200000 ,MAXALLIN=250000)
+ PARAMETER (MAXPTC=50000)
+ PARAMETER ( MAXFASI=1000, MAXBANDE=5)
+ PARAMETER (MAXFLAG=21 , MAXPAR=128 , MAXPARS=20 )
+ PARAMETER (MAXCELL=1000 )
+ PARAMETER (MAXTITLE=5)
+ PARAMETER ( MAXFILT=5 , MAXLAM=15 )
+C
+ COMMON /SURF/ NPNTMX,NPNT,NPNT1,NPNT2,NPNT3,
+ A N1I,N1F,N2I,N2F,N3I,N3F,
+ 1 X(MAXPT),Y(MAXPT),Z(MAXPT),
+ 2 G(MAXPT),GX(MAXPT),GY(MAXPT),GZ(MAXPT),
+ 3 FX(MAXPT),FY(MAXPT),FZ(MAXPT),
+ 4 TX(MAXPT),TY(MAXPT),TZ(MAXPT),
+ 5 T(MAXPT),A(MAXPT),
+ 6 FB(MAXPT),F(MAXPT,MAXBANDE),FRIFL(MAXPT),
+ 7 XLIMB(MAXPT),NUMOBJ(MAXPT),ALB(MAXPT),
+ 8 NUMCOARSE(MAXPT)
+C
+ DIMENSION NP123(3)
+ DIMENSION NIF(2,3)
+ EQUIVALENCE (NP123(1),NPNT1),(NIF(1,1),N1I)
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+C
+ COMMON /ALLIN/ NALLMX,NALL,ITERDONE,
+ 1 ISOUR(MAXALLIN),JRIC(MAXALLIN),
+ 2 TRANSFI(MAXALLIN),TRANSFJ(MAXALLIN),
+ 3 FSOURI(MAXALLIN),FSOURJ(MAXALLIN),
+ 4 FRICI(MAXALLIN),FRICJ(MAXALLIN),
+ 5 RIJ(MAXALLIN),
+ 6 RXIJ(MAXALLIN),RYIJ(MAXALLIN),RZIJ(MAXALLIN),
+ 7 COSGI(MAXALLIN),COSGJ(MAXALLIN),
+ 8 IAUS(MAXALLIN)
+C
+ COMMON /LIGHT/ NFASIMX,NFASI,NBANDE,
+ 1 FASE(MAXFASI),O(MAXBANDE,MAXFASI),
+ 2 C(3,MAXBANDE,MAXFASI),CT(MAXBANDE,MAXFASI),
+ 3 CBOL(3,MAXFASI),CBOLT(MAXFASI),
+ 4 CBOLR(3,MAXFASI),CBOLTR(MAXFASI),
+ 5 AREA(3),FBOLTOT(3),FBOL(3),FBOLREF(3)
+C
+ COMMON /VISUAL/ NCELL,CELLMAX,CELLMIN,DCELL,
+ 1 XCELL(MAXCELL),YCELL(MAXCELL)
+C
+ COMMON/TITLE/NTITLEMX,NTITLE,TITLE(MAXTITLE)
+ CHARACTER*80 TITLE
+C
+ COMMON /ROCHE/ AL1,AL2,AL3,XA1,XA2,XB1,XB2,YA,YB,
+ 1 OMEGAL1,OMEGAL2,OMEGAL3,
+ 2 RPOLE(3),XPOLE(3),YPOLE(3),ZPOLE(3),TPOLE(3),
+ 3 GXPOLE(3),GYPOLE(3),GZPOLE(3)
+C
+ COMMON /FILTRI/ NFILT,NLAM, NLAMBDA(MAXFILT),DELTALAM(MAXLAM),
+ 1 ALAMBDA(MAXLAM,MAXFILT),TRASMISS(MAXLAM,MAXFILT),
+ 2 COMPFRAC(MAXFILT,3)
+C
+ COMMON /PAR/ NFLAG,IFLAG(MAXFLAG),NPAR,PAR(MAXPAR)
+ DIMENSION PARS(MAXPARS,3)
+ EQUIVALENCE (PARS(1,1),PAR(1))
+ DIMENSION IPPAR1(MAXFLAG),IPPAR2(MAXFLAG)
+ CHARACTER*(*) NAMFLAG(MAXFLAG),DESCFLAG(MAXFLAG)
+ CHARACTER*(*) NAMPAR(MAXPAR),DESCPAR(MAXPAR)
+C
+ DIMENSION AUS(MAXPT)
+C
+ DIMENSION FASI01(MAXFASI)
+ CHARACTER*3 YESNO(-1:2)
+ DATA YESNO/'OLD','NO ','YES','NEW'/
+C primo: flag used to know if some line has been printed
+ LOGICAL PRIMO
+C ............................................................
+C
+C ............... phases from 0 to 1 for printing only
+ DO 10 I=1,NFASI
+ 10 FASI01(I)=FASE(I)/(3.14159*2.)
+C
+C ............... computing radii for each surf element if needed
+ IF(AMOUNT.GE.6..AND.(NT6.GT.0.OR.(NT.GT.0.AND.NT.LE.99) )) THEN
+ DO 20 I=1,NPNT
+ NUMERO=NUMOBJ(I)
+ X000=PARS(3,NUMERO)
+ Y000=PARS(4,NUMERO)
+ Z000=PARS(5,NUMERO)
+ AUS(I)=SQRT( (X(I)-X000)*(X(I)-X000)+
+ 1 (Y(I)-Y000)*(Y(I)-Y000)+
+ 2 (Z(I)-Z000)*(Z(I)-Z000) )
+ 20 CONTINUE
+ ENDIF
+C
+C ------------------------------ PRINTING ON UNIT 6
+ IF(NT6.GT.0) THEN
+C
+ IF(AMOUNT.LT.1.0) GO TO 600
+C ......................... general and input data
+ WRITE(6,5800)
+ 5800 FORMAT(1H1///,50X,29('-')/50X,' Ligth Curve Analysis Results'/
+ 1 50X,29('-')// )
+ WRITE(6,5799) (TITLE(J),J=1,NTITLE)
+ 5799 FORMAT(20X,A80)
+ WRITE(6,5801)
+ 5801 FORMAT(//' OPTION SUMMARY:')
+ WRITE(6,5802) (NAMFLAG(J),DESCFLAG(J),YESNO(IFLAG(J)),J=1,10)
+ 5802 FORMAT(1X,A,A,4X,A)
+C .................... Printing orbital parameters
+ WRITE(6,5813)
+ 5813 FORMAT(//' ORBITAL AND GRID PARAMETERS:')
+ K1=IPPAR1(11)
+ K2=IPPAR2(11)
+ WRITE(6,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ K1=IPPAR1(12)
+ K2=IPPAR2(12)
+ WRITE(6,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ K1=IPPAR1(13)
+ K2=IPPAR2(13)
+ WRITE(6,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ 5804 FORMAT(1X,A,A,4X,G15.7)
+C ..................... Printing reflection parameters
+ IF(IFLAG(9).GT.0) THEN
+ WRITE(6,5805)
+ 5805 FORMAT(/' PARAMETERS FOR REFLECTION COMPUTATION')
+ K1=IPPAR1(9)
+ K2=IPPAR2(9)
+ WRITE(6,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ ENDIF
+C .......................... printing colour band parameter
+ DO 30 I=4,7
+ IF(IFLAG(I).LE.0) GOTO 30
+ K1=IPPAR1(I)
+ K2=IPPAR2(I)
+ WRITE(6,5830) I-3,(NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ 5830 FORMAT(/' INPUT PARAMETERS FOR COLOUR BAND: ',I2/
+ 1 ((2X,A,2X,A,2X,G15.7)) )
+ 30 CONTINUE
+C .......................... printing print input parameters
+ K1=IPPAR1(20)
+ K2=IPPAR2(20)
+ WRITE(6,6015) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ 6015 FORMAT(//' PRINTED OUTPUT PARAMETERS:'/((2X,A,2X,A,2X,G15.7)))
+C
+C .......................... printing plot input parameters
+ IF(IFLAG(21).GT.0) THEN
+ K1=IPPAR1(21)
+ K2=IPPAR2(21)
+ WRITE(6,6017) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ 6017 FORMAT(//' PLOTTING OUTPUT PARAMETERS:'/((2X,A,2X,A,2X,G15.7)))
+ ENDIF
+C
+C .................. Printing stars and disk parameters
+ DO 60 I=1,3
+ IF(NP123(I).LE.0) GOTO 60
+ K1=(I-1)*MAXPARS+1
+ K2=I*MAXPARS
+ WRITE(6,6000) I,(NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ 6000 FORMAT(/' INPUT PARAMETERS FOR OBJECT NUMBER: ',I2/
+ 1 ((2X,A,2X,A,2X,G15.7)) )
+ 60 CONTINUE
+C .............. roche model parameters
+ IF(PARS(1,1).EQ.2..OR.PARS(1,2).EQ.2.) THEN
+ WRITE(6,6005) AL1,AL2,AL3,OMEGAL1,OMEGAL2,OMEGAL3,
+ 1 XA1,XA2,XB1,XB2,YA,YB
+ 6005 FORMAT(/' ROCHE MODEL DATA:'/
+ 1 ' L1,L2,L3 : lagrange points x position:',3G15.7/
+ 2 ' potential for surfaces at L1,L2,L3 :',3G15.7/
+ 3 ' XA1,XA2 : intersection star A-x axis:',2G15.7/
+ 4 ' XB1,XB2 : intersection star B-x axis:',2G15.7/
+ 5 ' approx. y dimension for critical lobe:',2G15.7)
+ PRIMO=.TRUE.
+ DO 62 I=1,3
+ IF(PARS(1,I).NE.2..OR.PARS(9,I).LE.0.0) GOTO 62
+C Following data is computed only for roche model and gravity dark.
+ IF(PRIMO) THEN
+ WRITE(6,6006)
+ 6006 FORMAT(//' Obj. ,pole: R',7X,'X',7X,7X,'Y',7X,7X,'Z',7X,
+ 1 7X,'GX',6X,7X,'GY',6X,7X,'GZ',6X,10X,'T')
+ PRIMO=.FALSE.
+ ENDIF
+ WRITE(6,6007) I,RPOLE(I),XPOLE(I),YPOLE(I),ZPOLE(I),
+ 2 GXPOLE(I),GYPOLE(I),GZPOLE(I),TPOLE(I)
+ 6007 FORMAT(1X,I2,1X,8G15.7)
+ 62 CONTINUE
+ ENDIF
+C ........................... run data
+ K1=INT(PAR(IPPAR1(9)+3))
+ WRITE(6,5803) NPNT,NPNTMX,NPNT1,NPNT2,NPNT3,
+ 1 NPNTC,NPNTMXC,K1,NPNT1C,NPNT2C,NPNT3C,
+ 2 NALL,MAXALLIN,ITERDONE
+ 5803 FORMAT(/' RUN PARAMETERS:'//
+ 1 10X,' Fine grid surface elements :',I10,10X,' Max allowed :',
+ 2 I10/48X,' For object 1:',I10/48X,' For object 2:',I10/
+ 3 48X,' For object 3:',I10//
+ 4 10X,' Coarse grid surface elements:',I10,10X,' Max allowed :',
+ 5 I10,' Coarsing factor:',I5/
+ 6 48X,' For object 1:',I10/48X,' For object 2:',I10/
+ 6 48X,' For object 3:',I10/
+ 7 10X,' Reflection paths:',I10,10X,' Max allowed :',I10,
+ 8 10X,' Iterations done:',I5)
+C
+ WRITE(6,6010)
+ 6010 FORMAT(///' BOLOMETRIC LUMINOSITY FOR EACH OBJECT :'/
+ 1 ' object number, total surf.area,',
+ 2 ' emitted lumin.,reflected lum.,'
+ 3 ' total lum., input lum. value, lum.norm.factor')
+ DO 61 I=1,3
+ IF(NP123(I).LE.0) GOTO 61
+ FBOLTR=FBOL(I)+FBOLREF(I)
+ WRITE(6,6011) I,AREA(I),FBOL(I),FBOLREF(I),FBOLTR,
+ 1 PARS(11,I),FBOLTOT(I)
+ 6011 FORMAT(I9,7X,6G15.7)
+ 61 CONTINUE
+C
+ WRITE(6,6012)
+ WRITE(6,6013) (I,(COMPFRAC(IB,I),IB=1,NFILT),I=1,3)
+ 6012 FORMAT(//' Objectc ',7X,'U-band lum.',5X,'B-band lum.',5X,
+ 1 'V-band lum.',5X,'R-band lum.',5X,'I-band lum.')
+ 6013 FORMAT((1X,I4,5X,5(1X,G15.7)))
+C
+ IF(AMOUNT.LT.2.0) GOTO 600
+C ...................... light received for each phase
+ WRITE(6,6020) NCELL,NCELL,DCELL,CELLMIN,CELLMAX
+ 6020 FORMAT(///' LIGHT CURVE PRODUCED:',8X,
+ 1 '(grid used:',I4,' *',I4,5X,'; step:',G12.5,2X,' from:',
+ 2 G12.5,5X,' to:',G12.5,' )'/5X,'Phase',5X,13X,
+ 3 5X,'bol.lum.',6X,'from star A',4X,'from star B',5X,'from disk')
+ WRITE(6,6021)
+ 1 (J,FASE(J),FASI01(J),CBOLT(J),CBOL(1,J),CBOL(2,J),CBOL(3,J),
+ 2 J=1,NFASI)
+ 6021 FORMAT((I4,G12.5,1X,G12.5,4(1X,G14.7)))
+C
+ IF(AMOUNT.LT.2.5) GOTO 600
+C .......................... reflected light for each fase
+ IF(IFLAG(9).GT.0) THEN
+ WRITE(6,6025)
+ 6025 FORMAT(///' REFLECTED LIGHT:'/5X,'Phase',5X,13X,
+ 3 4X,'refl.lum.',6X,'from star A',4X,'from star B',5X,'from disk')
+ WRITE(6,6021)
+ 1 (J,FASE(J),FASI01(J),CBOLTR(J),CBOLR(1,J),CBOLR(2,J),CBOLR(3,J),
+ 2 J=1,NFASI)
+C
+C .......................... light curve without reflection
+ WRITE(6,6026)
+ 6026 FORMAT(///' LIGHT CURVE WITHOUT REFLECTION:'/5X,'Phase',5X,13X,2X
+ 3,'unrefl.lum.',6X,'from star A',4X,'from star B',5X,'from disk')
+ WRITE(6,6021)
+ 1 (J,FASE(J),FASI01(J),CBOLT(J)-CBOLTR(J),
+ 2 CBOL(1,J)-CBOLR(1,J),CBOL(2,J)-CBOLR(2,J),
+ 3 CBOL(3,J)-CBOLR(3,J),J=1,NFASI)
+ ENDIF
+C
+ IF(AMOUNT.LT.3.0) GOTO 600
+C ...................... Color band flux
+ DO 65 IB=1,NBANDE
+ IF(IFLAG(3+IB).LE.0) GOTO 65
+ WRITE(6,6030) IB
+ 6030 FORMAT(///' LIGHT CURVE FOR COLOR ',I2,' :'/5X,'Phase',5X,13X,2X
+ 3,'luminosity ',6X,'from star A',4X,'from star B',5X,'from disk',
+ 4 4X,' observed ')
+ WRITE(6,6031)
+ 1 (J,FASE(J),FASI01(J),CT(IB,J),C(1,IB,J),C(2,IB,J),C(3,IB,J),
+ 2 O(IB,J),J=1,NFASI)
+ 6031 FORMAT((I4,G12.5,1X,G12.5,5(1X,G14.7) ))
+ 65 CONTINUE
+C
+ IF(AMOUNT.LT.5.0) GOTO 600
+C ..................... Coarse grid description of surfaces
+ IF(IFLAG(9).GT.0) THEN
+ WRITE(6,6045)
+ 6045 FORMAT(///' COARSE GRID SURFACE ELEMENTS DESCRIBING OBJECTS:'
+ a /' N. object',
+ 1 4X,'x',13X,'y',12X,'z',11X,'area',
+ 2 12X,'g',13X,'gx',11X,'gy',12X,'gz',11X,'num.fines'/)
+ WRITE(6,6046)
+ 1 (J,NUMOBJC(J),XC(J),YC(J),ZC(J),AC(J),
+ 2 GC(J),GXC(J),GYC(J),GZC(J),NUMFINI(J),J=1,NPNTC)
+ 6046 FORMAT((1X,2I3,8(1X,1PG13.6),1X,I9))
+ WRITE(6,6047)
+ 6047 FORMAT(///3X,'N. object',3X,' tot.bol.flux',
+ 1 5X,'reflected',8X,'limb. dark.',7X,'albedo',
+ 2 10X,' X-dim.',10X,' Y-dim.',10X,' Z-dim.')
+ WRITE(6,6048)
+ 1 (J,NUMOBJC(J),FBC(J),FRIFLC(J),XLIMBC(J),ALBC(J),
+ 2 DIMCX(J),DIMCY(J),DIMCZ(J),J=1,NPNTC)
+ 6048 FORMAT((1X,2I4,7(2X,G15.7)))
+ ENDIF
+C
+ IF(AMOUNT.LT.6.0) GOTO 600
+C ..................... Fine grid description of surfaces
+ WRITE(6,6050)
+ 6050 FORMAT(///' FINE GRID SURFACE ELEMENTS DESCRIBING OBJECTS:'
+ a /' N. object',
+ 1 3X,'x',13X,'y',12X,'z',11X,'area',3X,'coarse point'
+ 2 6X,'g',13X,'gx',13X,'gy',10X,'gz'/)
+ WRITE(6,6051)
+ 1 (J,NUMOBJ(J),X(J),Y(J),Z(J),A(J),NUMCOARSE(J),
+ 2 G(J),GX(J),GY(J),GZ(J),J=1,NPNT)
+ 6051 FORMAT((1X,I6,I2,4(1X,1PG13.6),1X,I6,4(1X,1PG13.6)))
+ WRITE(6,6052)
+ 6052 FORMAT(///3X,'N. object',3X,'temperature',3X,' tot.bol.flux',
+ 1 8X,'reflected',8X,'limb. dark.',7X,'albedo',10X,'radius')
+ WRITE(6,6053)
+ 1(J,NUMOBJ(J),T(J),FB(J),FRIFL(J),XLIMB(J),ALB(J),AUS(J),J=1,NPNT)
+ 6053 FORMAT((1X,I5,I3,6(2X,1PG15.7)))
+ WRITE(6,6054)
+ 6054 FORMAT(///5X,' ESTENSION FOR EACH SURFACE ELEMENT:'//
+ 1 3X,'N. object',8X,'Fx',13X,'Fy',17X,'Fz',
+ 2 13X,'Tx',16X,'Ty',15X,'Tz')
+ WRITE(6,6055)
+ 1 (J,NUMOBJ(J),FX(J),FY(J),FZ(J),TX(J),TY(J),TZ(J),J=1,NPNT)
+ 6055 FORMAT((1X,I5,I3,6(2X,1PG15.7)))
+C ......................
+C flux for each band : omissis
+C
+ IF(AMOUNT.LT.7) GOTO 600
+C .......................... Alignement for reflection
+ IF(IFLAG(9).LE.0) GOTO 600
+ WRITE(6,6070) ITERDONE,PAR(IPPAR1(9))
+ 6070 FORMAT(///' REFLECTION PATHS:'/
+ 1 40X,' iterations:',I5,5X,' albedo:',G15.7/
+ 2 ' path,sour-rec.coar.point,transf.func.,inv.trans,',
+ 3 ' d**2',7X,'cosgi',7X,'cosgj',2x,'last iter.source from I,J',
+ 4 ' last iter. received by I,J')
+ WRITE(6,6071) (J,ISOUR(J),JRIC(J),TRANSFI(J),TRANSFJ(J),RIJ(J),
+ 1COSGI(J),COSGJ(J),FSOURI(J),FSOURJ(J),FRICI(J),FRICJ(J),J=1,NALL)
+ 6071 FORMAT((1X,I5,2I6,4X,9(1X,G11.4)))
+C
+ ENDIF
+ 600 CONTINUE
+C
+C ------------------------------ PRINTING ON UNIT NT
+ IF(NT.GT.0.AND.NT.LE.99) THEN
+ IF(AMOUNT.LT.1.0) GO TO 700
+C ......................... general and input data
+ WRITE(NT,5800)
+ WRITE(NT,5799) (TITLE(J),J=1,NTITLE)
+ WRITE(NT,5801)
+ WRITE(NT,5802) (NAMFLAG(J),DESCFLAG(J),YESNO(IFLAG(J)),J=1,10)
+C Printing orbital parameters
+ WRITE(NT,5813)
+ K1=IPPAR1(11)
+ K2=IPPAR2(11)
+ WRITE(NT,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ K1=IPPAR1(12)
+ K2=IPPAR2(12)
+ WRITE(NT,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ K1=IPPAR1(13)
+ K2=IPPAR2(13)
+ WRITE(NT,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+C ........................Printing reflection parameters
+ IF(IFLAG(9).GT.0) THEN
+ WRITE(NT,5805)
+ K1=IPPAR1(9)
+ K2=IPPAR2(9)
+ WRITE(NT,5804) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ ENDIF
+C .......................... printing colour band parameter
+ DO 31 I=4,7
+ IF(IFLAG(I).LE.0) GOTO 31
+ K1=IPPAR1(I)
+ K2=IPPAR2(I)
+ WRITE(NT,5830) I-3,(NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ 31 CONTINUE
+C .......................... printing plot input parameters
+ IF(IFLAG(21).GT.0) THEN
+ K1=IPPAR1(21)
+ K2=IPPAR2(21)
+ WRITE(NT,6017) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+ ENDIF
+C .......................... printing print input parameters
+ K1=IPPAR1(20)
+ K2=IPPAR2(20)
+ WRITE(NT,6015) (NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+C
+C ........................Printing stars and disk parameters
+ DO 70 I=1,3
+ IF(NP123(I).LE.0) GOTO 70
+ K1=(I-1)*MAXPARS+1
+ K2=I*MAXPARS
+ WRITE(NT,6000) I,(NAMPAR(J),DESCPAR(J),PAR(J),J=K1,K2)
+C WRITE(NT,6000) I,(PARS(JJ,I),JJ=1,13)
+ 70 CONTINUE
+C ....................... roche model parameters
+ IF(PARS(1,1).EQ.2..OR.PARS(1,2).EQ.2.) THEN
+ WRITE(NT,6005) AL1,AL2,AL3,OMEGAL1,OMEGAL2,OMEGAL3,
+ 1 XA1,XA2,XB1,XB2,YA,YB
+ PRIMO=.TRUE.
+ DO 72 I=1,3
+ IF(PARS(1,I).NE.2..OR.PARS(9,I).LE.0.0) GOTO 72
+C Following data is computed only for roche model and gravity dark.
+ IF(PRIMO) THEN
+ WRITE(NT,6006)
+ PRIMO=.FALSE.
+ ENDIF
+ WRITE(NT,6007) I,RPOLE(I),XPOLE(I),YPOLE(I),ZPOLE(I),
+ 2 GXPOLE(I),GYPOLE(I),GZPOLE(I),TPOLE(I)
+ 72 CONTINUE
+ ENDIF
+C ........................... run data
+ K1=INT(PAR(IPPAR1(9)+3))
+ WRITE(NT,5803) NPNT,NPNTMX,NPNT1,NPNT2,NPNT3,
+ 1 NPNTC,NPNTMXC,K1,NPNT1C,NPNT2C,NPNT3C,
+ 2 NALL,MAXALLIN,ITERDONE
+C
+ WRITE(NT,6010)
+ DO 71 I=1,3
+ IF(NP123(I).LE.0) GOTO 71
+ FBOLTR=FBOL(I)+FBOLREF(I)
+ WRITE(NT,6011) I,AREA(I),FBOL(I),FBOLREF(I),FBOLTR,
+ 1 PARS(11,I),FBOLTOT(I)
+ 71 CONTINUE
+C
+ WRITE(NT,6012)
+ WRITE(NT,6013) (I,(COMPFRAC(IB,I),IB=1,NFILT),I=1,3)
+C
+ IF(AMOUNT.LT.2.0) GOTO 700
+C ...................... light received for each phase
+ WRITE(NT,6020) NCELL,NCELL,DCELL,CELLMIN,CELLMAX
+ WRITE(NT,6021)
+ 1 (J,FASE(J),FASI01(J),CBOLT(J),CBOL(1,J),CBOL(2,J),CBOL(3,J),
+ 2 J=1,NFASI)
+C
+ IF(AMOUNT.LT.2.5) GOTO 600
+C .......................... reflected light for each fase
+ IF(IFLAG(9).GT.0) THEN
+ WRITE(NT,6025)
+ WRITE(NT,6021)
+ 1 (J,FASE(J),FASI01(J),CBOLTR(J),CBOLR(1,J),CBOLR(2,J),CBOLR(3,J),
+ 2 J=1,NFASI)
+C .......................... light curve without reflection
+ WRITE(NT,6026)
+ WRITE(NT,6021)
+ 1 (J,FASE(J),FASI01(J),CBOLT(J)-CBOLTR(J),
+ 2 CBOL(1,J)-CBOLR(1,J),CBOL(2,J)-CBOLR(2,J),
+ 3 CBOL(3,J)-CBOLR(3,J),J=1,NFASI)
+ ENDIF
+C
+ IF(AMOUNT.LT.3.0) GOTO 700
+C ...................... Color band flux
+ DO 75 IB=1,NBANDE
+ IF(IFLAG(3+IB).LE.0) GOTO 75
+ WRITE(NT,6030) IB
+ WRITE(NT,6031)
+ 1 (J,FASE(J),FASI01(J),CT(IB,J),C(1,IB,J),C(2,IB,J),C(3,IB,J),
+ 2 O(IB,J),J=1,NFASI)
+ 75 CONTINUE
+C
+ IF(AMOUNT.LT.5.0) GOTO 700
+C ..................... Coarse grid description of surfaces
+ IF(IFLAG(9).GT.0) THEN
+ WRITE(NT,6045)
+ WRITE(NT,6046)
+ 1 (J,NUMOBJC(J),XC(J),YC(J),ZC(J),AC(J),
+ 2 GC(J),GXC(J),GYC(J),GZC(J),NUMFINI(J),J=1,NPNTC)
+ WRITE(NT,6047)
+ WRITE(NT,6048)
+ 1 (J,NUMOBJC(J),FBC(J),FRIFLC(J),XLIMBC(J),ALBC(J),
+ 2 DIMCX(J),DIMCY(J),DIMCZ(J),J=1,NPNTC)
+ ENDIF
+C
+ IF(AMOUNT.LT.6.0) GOTO 700
+C ..................... Fine grid Description of surfaces
+ WRITE(NT,6050)
+ WRITE(NT,6051)
+ 1 (J,NUMOBJ(J),X(J),Y(J),Z(J),A(J),NUMCOARSE(J),
+ 2 G(J),GX(J),GY(J),GZ(J),J=1,NPNT)
+ WRITE(NT,6052)
+ WRITE(NT,6053)
+ 1(J,NUMOBJ(J),T(J),FB(J),FRIFL(J),XLIMB(J),ALB(J),AUS(J),J=1,NPNT)
+ WRITE(NT,6054)
+ WRITE(NT,6055)
+ 1 (J,NUMOBJ(J),FX(J),FY(J),FZ(J),TX(J),TY(J),TZ(J),J=1,NPNT)
+C ......................
+C flux for each band : omissis
+C
+ IF(AMOUNT.LT.7) GOTO 700
+C .......................... Alignement for reflection
+ IF(IFLAG(9).LE.0) GOTO 700
+ WRITE(NT,6070) ITERDONE,PAR(IPPAR1(9))
+ WRITE(NT,6071) (J,ISOUR(J),JRIC(J),TRANSFI(J),TRANSFJ(J),RIJ(J),
+ 1COSGI(J),COSGJ(J),FSOURI(J),FSOURJ(J),FRICI(J),FRICJ(J),J=1,NALL)
+C
+ ENDIF
+ 700 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE REFLECT(NITER,PKCONV,R4PREC,ALBEDO,
+ 1 NPRINTR,NFILE,FORDER)
+C ------------------------------------------------
+C Computes the light reflected by the 3 objects,
+C accounting for the geometry of stars and disk.
+C Reflection of increasing orders is computed,
+C until the order NITER or if the convergence
+C criterion specified by PKCONV is met.
+C ALBEDO is an input parameter specifing the
+C efficiency of the reflection process
+C
+C R4PREC: sin,cos < R4prec are assumed 0 to shorten computations
+C and avoid rounding error in sin,cos whithout using double precision.
+C ------------------------------------------------
+ PARAMETER (MAXPT=200000 ,MAXALLIN=250000, MAXBANDE=5)
+ PARAMETER (MAXPTC=50000)
+C
+ COMMON /SURFC/ NPNTMXC,NPNTC,NPNT1C,NPNT2C,NPNT3C,
+ A N1IC,N1FC,N2IC,N2FC,N3IC,N3FC,
+ 1 XC(MAXPTC),YC(MAXPTC),ZC(MAXPTC),
+ 2 GC(MAXPTC),GXC(MAXPTC),GYC(MAXPTC),GZC(MAXPTC),
+ 5 AC(MAXPTC),
+ 6 DIMCX(MAXPTC),DIMCY(MAXPTC),DIMCZ(MAXPTC),
+ 6 FBC(MAXPTC),FRIFLC(MAXPTC),
+ 7 XLIMBC(MAXPTC),NUMOBJC(MAXPTC),ALBC(MAXPTC),
+ 8 NUMFINI(MAXPTC)
+ DIMENSION NP123C(3)
+ DIMENSION NIFC(2,3)
+ EQUIVALENCE (NP123C(1),NPNT1C),(NIFC(1,1),N1IC)
+C
+ COMMON /ALLIN/ NALLMX,NALL,ITERDONE,
+ 1 ISOUR(MAXALLIN),JRIC(MAXALLIN),
+ 2 TRANSFI(MAXALLIN),TRANSFJ(MAXALLIN),
+ 3 FSOURI(MAXALLIN),FSOURJ(MAXALLIN),
+ 4 FRICI(MAXALLIN),FRICJ(MAXALLIN),
+ 5 RIJ(MAXALLIN),
+ 6 RXIJ(MAXALLIN),RYIJ(MAXALLIN),RZIJ(MAXALLIN),
+ 7 COSGI(MAXALLIN),COSGJ(MAXALLIN),
+ 8 IAUS(MAXALLIN)
+C
+ DIMENSION FORDER(MAXPTC)
+C pi4=1/(4*pi)
+ PARAMETER PI4=0.07957747
+C ...........................................................
+C Zero alignements at the beginning
+ NALL=0
+C Find allignements (possible reflection paths)
+C between objects 1 and 2 ( star A and star B)
+C An alignement list is built in common /allin/
+ CALL FINDALL(1,2,R4PREC)
+C First and last alignement betwenn object 1 and 2 (A and B)
+ NALLAB1=1
+ NALLAB2=NALL
+C Eliminates allignements between stars A and B
+C which are intercepted by the disk (object 3, or C)
+C Common /allin/ is updated.
+ CALL ALLSEL(3,NALLAB1,NALLAB2)
+ NALLAB2=NALL
+C Find allignements between objects 1 and 3
+C ( star A and disk ), the alignement list
+C in common /allin/ is updated
+ CALL FINDALL(1,3,R4PREC)
+ NALLAC1=NALLAB2+1
+ NALLAC2=NALL
+C Eliminates allignement between A and C intercepted by B
+ CALL ALLSEL(2,NALLAC1,NALLAC2)
+ NALLAC2=NALL
+C Find allignements between objects 2 and 3
+C ( star B and disk ), the alignement list
+C in common /allin/ is updated
+ CALL FINDALL(2,3,R4PREC)
+ NALLBC1=NALLAC2+1
+ NALLBC2=NALL
+C Eliminates the alignements between B and C intercepted by A
+ CALL ALLSEL(1,NALLBC1,NALLBC2)
+ NALLBC2=NALL
+C Completing the data in the alignement list.
+C computing FSOUR, initial alignement source
+ CALL CALCALL(ALBEDO)
+C Zero reflection source for each surface element
+ CALL NULL(NPNTC,FRIFLC)
+C Initial value for convergency checks
+ FMAX=0.0
+ DO 40 K=1,NALL
+ FMAX=FMAX+FSOURI(K)+FSOURJ(K)
+ 40 CONTINUE
+C ------------------------- LOOP ON REFLECTION ORDERS
+ DO 50 IT=1,NITER
+C
+C Computes the light transmitted along each allignement
+ DO 60 K=1,NALL
+ FRICJ(K)=FSOURI(K)*TRANSFI(K)
+ 60 FRICI(K)=FSOURJ(K)*TRANSFJ(K)
+C Zeroes the light falling on on each coarse surf.elem. for this order
+ DO 55 K=1,NPNTC
+ 55 FORDER(K)=0.0
+C Computes the light received by each coarse surface element
+C At this reflection order (forder) and the total FRIFLC
+C Summing on all path from I to J and inverse (J to I)
+ DO 70 K=1,NALL
+ FORDER(JRIC(K)) =FORDER(JRIC(K)) +FRICJ(K)
+ FORDER(ISOUR(K))=FORDER(ISOUR(K))+FRICI(K)
+ FRIFLC(JRIC(K)) =FRIFLC(JRIC(K)) +FRICJ(K)
+ 70 FRIFLC(ISOUR(K))=FRIFLC(ISOUR(K))+FRICI(K)
+C
+C Convergence check
+ FMAX1=0.0
+ DO 90 K=1,NALL
+ FMAX1=FMAX1+FRICI(K)+FRICJ(K)
+ 90 CONTINUE
+C
+C Prints reflection interation data
+ IF(NPRINTR.GT.0.AND.NFILE.GT.0.AND.NFILE.LE.99) THEN
+ WRITE(NFILE,8000) IT,NITER,PKCONV,FMAX,FMAX1
+ 8000 FORMAT(//' Reflection hystory: iteration: ',I3,' max iter:',
+ 1 I3,' conv.check:',G12.4,' initial,this iter refl.flux:',
+ 2 2G12.4/' path,I:sour.,J:rec.,',
+ 3 'flux from I, transf.func.I=>J,flux to J,'
+ 4 'flux from J, transf.func.J=>I,flux to I,'
+ 5 'Sum I:I=>J this iter.,all iter.')
+ DO 95 K=1,NALL
+ I=ISOUR(K)
+ J=JRIC(K)
+ WRITE(NFILE,8001) K,I,J,FSOURI(K),TRANSFI(K),FRICJ(K),
+ 1 FSOURJ(K),TRANSFJ(K),FRICJ(K),
+ 2 FORDER(K),FRIFLC(K)
+ 8001 FORMAT((1X,I5,2I6,8(1X,1P,G13.6)))
+ 95 CONTINUE
+ ENDIF
+C Convergency check continues:
+ IF(FMAX.NE.0.0) THEN
+ FFMAX=FMAX1/FMAX
+ ELSE
+ FFMAX=0.0
+ ENDIF
+C
+ IF(FFMAX.GT.1.)THEN
+ WRITE(6,1000) IT
+ IF(NFILE.GT.0.AND.NFILE.LE.99) WRITE(NFILE,1000) IT
+ 1000 FORMAT(' WARNING! Iteration:',I3,' REFLECTION DIVERGING')
+ ELSE IF(FFMAX.LT.PKCONV) THEN
+ ITERDONE=IT
+ GOTO 500
+ ENDIF
+C Recomputes sources for the next order of reflection
+ IF(IT.EQ.NITER) GOTO 50
+ DO 80 K=1,NALL
+ FSOURI(K)=FORDER(ISOUR(K))
+ 80 FSOURJ(K)=FORDER(JRIC(K))
+C
+ 50 CONTINUE
+ ITERDONE=NITER
+ 500 RETURN
+ END
+C
+ SUBROUTINE ROCHES(PRINT6,PRINTFILE,NFILE,NTH,NTHF,NCORD,NPUNT,
+ 1 NPUNTC, XA,XB,XSHIFT,OMEGA,Q,RIS,PRECISION,
+ 2 NPNTMX,X,Y,Z,G,GX,GY,GZ,FX,FY,FZ,TX,TY,TZ,AREA,A,NUMCOARSE,
+ 3 SINFI,COSFI,FUNCTION)
+C ----------------------------------------------------------------
+C Called for a Roche lobe star, fills the surface element list
+C with the surface elements of a roche lobe of given Omega.
+C The main star (object A) is build in (0,0,0) the star B in
+C (1,0,0) the star is then shifted on the x axis by XSHIFT.
+C RIS=0. for the main star A, otherwise 1.
+C XA is the maximum X value, XB the minimum one for this star;
+C for star A XA is near the L1 point, while for star B XB is near L1
+C FUNCTION is rochesf1 for main star A, rochesf2 for star B,
+C passed to RTNEWT, to compute the Radius value for a given angle
+C in the roche model. Q is the mass ratio;
+C SINFI,COSFI are working spaces.
+C
+C This routine has been structured in the same way as Sfera2;
+C but, due to a lack in equatorial simmetry for the Roche model,
+C the north part and the sud part are not treated togheter.
+C The arch lenght has been approximated by the sphere value:
+C NCORD=1: Theta and phi arch extension for each surface element
+C NCORD=2: Theta and phi cord extension of each surface element
+C NCORD=3: Theta and phi tangent segment of each surface element
+C T, and F tangent vectors are obtained in an approximate way
+C from the spherical ones; to shorten calculations
+C ----------------------------------------------------------------
+ COMMON /ROCHE/ AL1,AL2,AL3,XA1,XA2,XB1,XB2,YA,YB,
+ 1 OMEGAL1,OMEGAL2,OMEGAL3,
+ 2 RPOLE(3),XPOLE(3),YPOLE(3),ZPOLE(3),TPOLE(3),
+ 3 GXPOLE(3),GYPOLE(3),GZPOLE(3)
+ DIMENSION X(NPNTMX),Y(NPNTMX),Z(NPNTMX)
+ DIMENSION G(NPNTMX),GX(NPNTMX),GY(NPNTMX),GZ(NPNTMX)
+ DIMENSION FX(NPNTMX),FY(NPNTMX),FZ(NPNTMX)
+ DIMENSION TX(NPNTMX),TY(NPNTMX),TZ(NPNTMX)
+ DIMENSION A(NPNTMX)
+ DIMENSION NUMCOARSE(NPNTMX)
+ DIMENSION SINFI(NPNTMX),COSFI(NPNTMX)
+ LOGICAL PRINTFILE,PRINT6
+ EXTERNAL FUNCTION
+ DATA PI /3.1415926/
+C DATA PI /3.141592653589793/
+C
+C Newton-Cotes Error Counter
+ NERROR=0
+C Delta theta fine, number of phi fines (same as sfera2):
+ DTH=PI/(NTH-1)
+ DTH5=0.5*DTH
+C at equator: DFI=(2.*PI)/NFIEQ;DFI5=0.5*DFI then decreases as sin theta
+ NFIEQ=2*(NTH-1)
+C Theta coarse arc interval and number of coarse theta points
+ DTHC=DTH*NTHF
+ DTHC5=DTHC*0.5
+ NTHC=(NTH+(NTHF-1))/NTHF
+ IF(NTHC*NTHF.NE.NTH+NTHF-1) THEN
+ WRITE(6,900) NTH,NTHF
+ IF(PRINTFILE) WRITE(NFILE,900) NTH,NTHF
+ 900 FORMAT(/,' ERRROR! The number of theta points:',I6,
+ 1 ' is not consistent with the coarsing factor:',I6)
+ ENDIF
+C RMAX: Approximate max radius for cord calculation(for F and T):
+C NOW THE R VALUE FOR EACH THE THETA PHI IS USED INSTEAD OF RMAX!
+C RMAX=MAX(ABS(XA-RIS),ABS(XB-RIS))
+ IF(NCORD.LE.1) THEN
+C RTCORDA=DTH5*RMAX
+ RTCORDA=DTH5
+ ELSE IF(NCORD.EQ.2) THEN
+C RTCORDA=RMAX*SIN(DTH5)
+ RTCORDA=SIN(DTH5)
+ ELSE
+C RTCORDA=RMAX*TAN(DTH5)
+ RTCORDA=TAN(DTH5)
+ ENDIF
+C
+ NPUNT=0
+ ICT=NPUNTC
+ LOST=0
+ TH=0.0
+C Limits for Newton Cotes R-Search (star radius)
+ IF(RIS.EQ.0.0) THEN
+C Main star A
+ RLIM1=PRECISION
+ RLIM2=AL1
+ ELSE
+C Minor star B (this work if the max radius of B is in L1)
+ RLIM1=PRECISION
+ RLIM2=1.-AL1
+ ENDIF
+C ---------------
+C ................................ North polar cap
+C ---------------
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 500
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+C the first r value is the intersection of the lobe with x axis
+ R=XA-RIS
+ X1=XA
+ Y1=0.
+ Z1=0.
+C normal to the surf. in x1,y1,z1
+ IF(X1.NE.AL1.AND.X1.NE.AL2) THEN
+ CALL GRADOMEGA(X1,Y1,Z1,GX(NPUNT),GY(NPUNT),GZ(NPUNT),G(NPUNT),Q)
+ ELSE
+C Singular gradient in lagrangian point
+ X1MOD=X1-10.*PRECISION
+ CALL GRADOMEGA(X1MOD,Y1,Z1,GX(NPUNT),GY(NPUNT),GZ(NPUNT),
+ 1 G(NPUNT),Q)
+ GX(NPUNT)=1.
+ GY(NPUNT)=0.
+ GZ(NPUNT)=0.
+ ENDIF
+C cosine of angle between (x,y,z) and normal vector
+ COSB=1.
+ A(NPUNT)=(2.*PI)*(1.-COS(DTH5))*R*(R/COSB)
+ AREA=A(NPUNT)
+ X(NPUNT)=X1
+ Y(NPUNT)=Y1
+ Z(NPUNT)=Z1
+ FY(NPUNT)=0.
+ FZ(NPUNT)=RTCORDA*(R/COSB)
+ FX(NPUNT)=0.
+ TY(NPUNT)=RTCORDA*(R/COSB)
+ TZ(NPUNT)=0.
+ TX(NPUNT)=0.
+C the T and F component along G are subtracted, making T and F
+C perpendicular to G
+ CALL ROSUB(FX(NPUNT),FY(NPUNT),FZ(NPUNT),TX(NPUNT),TY(NPUNT),
+ 1 TZ(NPUNT),GX(NPUNT),GY(NPUNT),GZ(NPUNT),COSB)
+ NUMCOARSE(NPUNT)=ICT
+C ------------------------------
+C ...........................coarse cap: other theta points
+C ------------------------------
+ COSTH2=COS(DTH5)
+ NTHF2=NTHF/2+1
+ DO 5 IT=2,NTHF2
+ TH=DTH*(IT-1)
+ SINTH=SIN(TH)
+ COSTH=COS(TH)
+ COSTH1=COSTH2
+ COSTH2=COS(TH+DTH5)
+C .... Number of phi fines for this theta coarse
+ NFI=INT(NFIEQ*SINTH)
+ IF(NFI.LT.4) NFI=4
+ DFI=2.*PI/NFI
+ DFI5=DFI*0.5
+C .... F tangent Interval extension for this theta value
+C really i should use R instead of RMAX and put the following into loop 5
+ IF(NCORD.LE.1) THEN
+C RFCORDA=DFI5*RMAX
+ RFCORDA=DFI5
+ ELSE IF(NCORD.EQ.2) THEN
+C RFCORDA=RMAX*SIN(TH+DTH5)*SIN(DFI5)
+ RFCORDA=SIN(TH+DTH5)*SIN(DFI5)
+ ELSE
+C RFCORDA=RMAX*SIN(TH+DTH5)*TAN(DFI5)/COS(DTH5)
+ RFCORDA=SIN(TH+DTH5)*TAN(DFI5)/COS(DTH5)
+ ENDIF
+C
+ DO 5 IF=1,NFI
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 5
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ FI=DFI * ( IF-1-INT(NTHF/2) )
+ COSFII=COS(FI)
+ SINFII=SIN(FI)
+ SINCOS=SINTH*COSFII
+ SINSIN=SINTH*SINFII
+C find r for this theta,fi; preceeding R is given as a guess value
+C limits for r search are 0 and L1 (to avoid an r in the other object)
+ R=RTNEWT(FUNCTION,RLIM1,RLIM2,PRECISION,R,Q,OMEGA,COSTH,SINSIN,
+ 1 NERROR)
+ X1=R*COSTH+RIS
+ Y1=R*SINCOS
+ Z1=R*SINSIN
+ CALL GRADOMEGA(X1,Y1,Z1,GX(NPUNT),GY(NPUNT),GZ(NPUNT),G(NPUNT),Q)
+C for star B remember that cosb is the angle G-star center
+C COSB=((X1-RIS)*GX(NPUNT)+Y1*GY(NPUNT)+Z1*GZ(NPUNT))/R
+ COSB=(COSTH*GX(NPUNT)+SINCOS*GY(NPUNT)+SINSIN*GZ(NPUNT))
+ X(NPUNT)=X1
+ Y(NPUNT)=Y1
+ Z(NPUNT)=Z1
+ FY(NPUNT)=-SINFII*RFCORDA*(R/COSB)
+ FZ(NPUNT)= COSFII*RFCORDA*(R/COSB)
+ FX(NPUNT)=0.
+ TY(NPUNT)= (COSTH*RTCORDA)*COSFII*(R/COSB)
+ TZ(NPUNT)= (COSTH*RTCORDA)*SINFII*(R/COSB)
+ TX(NPUNT)=-(SINTH*RTCORDA)*(R/COSB)
+ CALL ROSUB(FX(NPUNT),FY(NPUNT),FZ(NPUNT),TX(NPUNT),TY(NPUNT),
+ 1 TZ(NPUNT),GX(NPUNT),GY(NPUNT),GZ(NPUNT),COSB)
+ DFIRR=DFI*R*R
+ AREATH=DFIRR*(COSTH1-COSTH2)
+ A(NPUNT)=AREATH/COSB
+ AREA=AREA+A(NPUNT)
+ NUMCOARSE(NPUNT)=ICT
+ 5 CONTINUE
+C ------------------
+C .................................... other theta points
+C ------------------
+C
+C ------------------------------ Loop on theta values
+C --------------------
+ NTH1=(NTHF/2)+2
+ NTH2=NTH-(NTHF/2)-1
+ TH1=(NTH1-1)*DTH
+ COSTH2=COS(TH1-DTH5)
+C num.of phi coarse for this theta to increment ICT for first step below
+ NFIC=1
+C
+ DO 20 IT=NTH1,NTH2
+ TH=(IT-1)*DTH
+ SINTH=SIN(TH)
+ COSTH=COS(TH)
+ COSTH1=COSTH2
+ COSTH2=COS(TH+DTH5)
+C ................ changing the theta coarse point each nthf fine points.
+ IF(MOD(IT+NTHF/2-1,NTHF).EQ.0) THEN
+ ICT=ICT+NFIC
+C Number of phi coarse for this nthf interval
+ THMEDIO=DTH*(IT-1+NTHF/2)
+ NFIC=INT(2.*PI*SIN(THMEDIO)/DTHC)
+ IF(NFIC.LT.4) NFIC=4
+C Number of phi fine for this coarse theta range
+ NFI=NFIC*NTHF
+ DFI=(2.*PI)/NFI
+ DFI5=DFI*0.5
+C Cos and sin fi values for this theta range computed once
+ DO 25 IF=1,NFI
+ FI=DFI * ( IF-1-INT(NTHF/2) )
+ COSFI(IF)=COS(FI)
+ SINFI(IF)=SIN(FI)
+ 25 CONTINUE
+ ENDIF
+C ................................
+C The first coarse num.in phi loop for this th. coar.is ICT
+ ICF=ICT
+C
+C really i should use R instead of RMAX and put the following into loop 20
+ IF(NCORD.LE.1) THEN
+C RFCORDA=DFI5*RMAX
+ RFCORDA=DFI5
+ ELSE
+C max diameter of sphere between TH+DTH5 and TH-DTH5: it is
+C +DTH5 over and -DTH5 below equator
+ IF(TH.GT.PI/2.) THEN
+ SINMAX=SIN(TH-DTH5)
+ ELSE
+ SINMAX=SIN(TH+DTH5)
+ ENDIF
+ IF(NCORD.EQ.2) THEN
+ RFCORDA=SINMAX*SIN(DFI5)
+ ELSE
+ RFCORDA=SINMAX*TAN(DFI5)/COS(DTH5)
+ ENDIF
+ ENDIF
+C ------------------------------------ Loop on phi values
+C --------------------
+ DO 20 IF=1,NFI
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 20
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ SINCOS=SINTH*COSFI(IF)
+ SINSIN=SINTH*SINFI(IF)
+ R=RTNEWT(FUNCTION,RLIM1,RLIM2,PRECISION,R,Q,OMEGA,COSTH,SINSIN,
+ 1 NERROR)
+ X1=R*COSTH+RIS
+ Y1=R*SINCOS
+ Z1=R*SINSIN
+ CALL GRADOMEGA(X1,Y1,Z1,GX(NPUNT),GY(NPUNT),GZ(NPUNT),G(NPUNT),Q)
+ COSB=(COSTH*GX(NPUNT)+SINCOS*GY(NPUNT)+SINSIN*GZ(NPUNT))
+C COSB=((X1-RIS)*GX(NPUNT)+Y1*GY(NPUNT)+Z1*GZ(NPUNT))/R
+ X(NPUNT)=X1
+ Y(NPUNT)=Y1
+ Z(NPUNT)=Z1
+ FY(NPUNT)=-SINFI(IF)*RFCORDA*(R/COSB)
+ FZ(NPUNT)= COSFI(IF)*RFCORDA*(R/COSB)
+ FX(NPUNT)=0.
+ TY(NPUNT)= (COSTH*RTCORDA)*COSFI(IF)*(R/COSB)
+ TZ(NPUNT)= (COSTH*RTCORDA)*SINFI(IF)*(R/COSB)
+ TX(NPUNT)=-(SINTH*RTCORDA)*(R/COSB)
+ CALL ROSUB(FX(NPUNT),FY(NPUNT),FZ(NPUNT),TX(NPUNT),TY(NPUNT),
+ 1 TZ(NPUNT),GX(NPUNT),GY(NPUNT),GZ(NPUNT),COSB)
+ DFIRR=DFI*R*R
+ AREATH=DFIRR*(COSTH1-COSTH2)
+ A(NPUNT)=AREATH/COSB
+ AREA=AREA+A(NPUNT)
+ NUMCOARSE(NPUNT)=ICF
+C AT each NTHF phi points change the coarse phi point
+C Including the remainder of the division IF/NTHF in the last coarse
+ IF(MOD(IF,NTHF).EQ.0) ICF=ICF+1
+C
+ 20 CONTINUE
+C
+C Next ICT to fill:
+ ICT=ICT+NFIC
+C ---------------
+C .................................... South pole
+C ---------------
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 500
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ R=RIS-XB
+ X1=XB
+ Y1=0.
+ Z1=0.
+C normal to the surf. in x1,y1,z1
+ IF(X1.NE.AL1.AND.X1.NE.AL3) THEN
+ CALL GRADOMEGA(X1,Y1,Z1,GX(NPUNT),GY(NPUNT),GZ(NPUNT),G(NPUNT),Q)
+ ELSE
+C Singular gradient in lagrangian point
+ X1MOD=X1+10.*PRECISION
+ CALL GRADOMEGA(X1MOD,Y1,Z1,GX(NPUNT),GY(NPUNT),GZ(NPUNT),
+ 1 G(NPUNT),Q)
+ GX(NPUNT)=-1.
+ GY(NPUNT)=0.
+ GZ(NPUNT)=0.
+ ENDIF
+C cosine of angle between (x,y,z) and normal vector
+ COSB=1.
+ A(NPUNT)=(2.*PI)*(1.-COS(DTH5))*R*R/COSB
+ AREA=AREA+A(NPUNT)
+ X(NPUNT)=X1
+ Y(NPUNT)=Y1
+ Z(NPUNT)=Z1
+ FY(NPUNT)=0.
+ FZ(NPUNT)=-RTCORDA*(R/COSB)
+ FX(NPUNT)=0.
+ TY(NPUNT)=RTCORDA*(R/COSB)
+ TZ(NPUNT)=0.
+ TX(NPUNT)=0.
+C the T and F component along G are subtracted, making T and F
+C perpendicular to G
+ CALL ROSUB(FX(NPUNT),FY(NPUNT),FZ(NPUNT),TX(NPUNT),TY(NPUNT),
+ 1 TZ(NPUNT),GX(NPUNT),GY(NPUNT),GZ(NPUNT),COSB)
+ NUMCOARSE(NPUNT)=ICT
+C ----------------------------------
+C ...........................sud coarse cap: other theta points
+C ----------------------------------
+ NTH1=NTH2+1
+ NTH2=NTH-1
+ TH1=(NTH1-1.)*DTH
+ COSTH2=COS(TH1-DTH5)
+ DO 45 IT=NTH1,NTH2
+ TH=DTH*(IT-1)
+ SINTH=SIN(TH)
+ COSTH=COS(TH)
+ COSTH1=COSTH2
+ COSTH2=COS(TH+DTH5)
+C .... Number of phi fines for this theta coarse
+ NFI=INT(NFIEQ*SINTH)
+ IF(NFI.LT.4) NFI=4
+ DFI=2.*PI/NFI
+ DFI5=DFI*0.5
+C .... F tangent Interval extension for this theta value
+C really i should use R instead of RMAX and put the following into loop 45
+ IF(NCORD.LE.1) THEN
+C RFCORDA=DFI5*R
+ RFCORDA=DFI5
+ ELSE IF(NCORD.EQ.2) THEN
+C RFCORDA=R*SIN(TH-DTH5)*SIN(DFI5)
+ RFCORDA=SIN(TH-DTH5)*SIN(DFI5)
+ ELSE
+C RFCORDA=R*SIN(TH-DTH5)*TAN(DFI5)/COS(DTH5)
+ RFCORDA=SIN(TH-DTH5)*TAN(DFI5)/COS(DTH5)
+ ENDIF
+C
+ DO 45 IF=1,NFI
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 5
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ FI=DFI * ( IF-1-INT(NTHF/2) )
+ COSFII=COS(FI)
+ SINFII=SIN(FI)
+ SINCOS=SINTH*COSFII
+ SINSIN=SINTH*SINFII
+C find r for this theta,fi; preceeding R is given as a guess value
+C limits for r search are 0 and L1 (to avoid an r in the other object)
+ R=RTNEWT(FUNCTION,RLIM1,RLIM2,PRECISION,R,Q,OMEGA,COSTH,SINSIN,
+ 1 NERROR)
+ X1=R*COSTH+RIS
+ Y1=R*SINCOS
+ Z1=R*SINSIN
+ CALL GRADOMEGA(X1,Y1,Z1,GX(NPUNT),GY(NPUNT),GZ(NPUNT),G(NPUNT),Q)
+ COSB=(COSTH*GX(NPUNT)+SINCOS*GY(NPUNT)+SINSIN*GZ(NPUNT))
+C COSB=((X1-RIS)*GX(NPUNT)+Y1*GY(NPUNT)+Z1*GZ(NPUNT))/R
+ X(NPUNT)=X1
+ Y(NPUNT)=Y1
+ Z(NPUNT)=Z1
+ FY(NPUNT)=-SINFII*RFCORDA*(R/COSB)
+ FZ(NPUNT)= COSFII*RFCORDA*(R/COSB)
+ FX(NPUNT)=0.
+ TY(NPUNT)= (COSTH*RTCORDA)*COSFII*(R/COSB)
+ TZ(NPUNT)= (COSTH*RTCORDA)*SINFII*(R/COSB)
+ TX(NPUNT)=-(SINTH*RTCORDA)*(R/COSB)
+ CALL ROSUB(FX(NPUNT),FY(NPUNT),FZ(NPUNT),TX(NPUNT),TY(NPUNT),
+ 1 TZ(NPUNT),GX(NPUNT),GY(NPUNT),GZ(NPUNT),COSB)
+ DFIRR=DFI*R*R
+ AREATH=DFIRR*(COSTH1-COSTH2)
+ A(NPUNT)=AREATH/COSB
+ AREA=AREA+A(NPUNT)
+ NUMCOARSE(NPUNT)=ICT
+ 45 CONTINUE
+C
+ NPUNTC=ICT+1
+C
+C Shifts star , if requested
+ IF(XSHIFT.NE.0.0) THEN
+ DO 50 J=1,NPUNT
+ 50 X(J)=X(J)+XSHIFT
+ ENDIF
+C
+ 500 IF(LOST.GT.0) THEN
+ WRITE(6,1000) NPUNT,LOST
+ IF(PRINTFILE) WRITE(NFILE,1000) NPUNT,LOST
+ 1000 FORMAT(' ERROR! OBJECT NOT PROPERLY DRAWN! '/
+ 1 ' space available for only:',I6,' surface points.',
+ 2 ' lost points:',I6/
+ 3 ' REDUCE THE NTHETA PARAMETER OR GO INTO THE FORTRAN LIST TO',
+ 4 ' INCREASE MAXPT')
+ ENDIF
+C
+ IF(NERROR.GT.0) THEN
+ WRITE(6,1100) NERROR
+ IF(PRINTFILE) WRITE(NFILE,1100) NERROR
+ 1100 FORMAT(/' ERROR : In roche lobe drawing: r not found for',
+ 1 I5,' surf.elements'/)
+ ENDIF
+ RETURN
+ END
+C
+ SUBROUTINE ROCHESF1(R,F,DF,Q,OMEGA,A,B)
+C ---------------------------------------------------------------
+C Used by Roches : lobe function and derivative for main object A
+C ---------------------------------------------------------------
+ DUM=1.+R*R-2.*R*A
+ SQR3=DUM**(-3./2.)
+ SQR=1./SQRT(DUM)
+C F=1./R-OMEGA+Q*SQR+0.5*(1.+Q)*R*R*(1.-B*B)-A*Q*R+0.5*Q*Q/(1.+Q)old def.
+C F=1./R-OMEGA+Q*SQR+0.5*(1.+Q)*R*R*(1.-B*B)-A*Q*R
+ F=1./R -OMEGA+ 0.5*(1.+Q)*(R*R)*(1.-B*B)+ Q*(SQR-A*R)
+ DF=(1.+Q)*R*(1.-B*B) -1./(R*R) -Q*( (R-A)*SQR3 + A)
+ RETURN
+ END
+C
+ SUBROUTINE ROCHESF2(R,F,DF,Q,OMEGA,A,B)
+C ---------------------------------------------------------------
+C Used by Roches : lobe function and derivative for minor object B
+C ---------------------------------------------------------------
+ DUM=1.+R*R+2.*R*A
+ SQR3=DUM**(-3./2.)
+ SQR=1./SQRT(DUM)
+C F=SQR-OMEGA+Q/R+R*A+0.5*(1.+Q)* R*R*(A*A+G*G) +0.5/(1.+Q) old omega def.
+ F=SQR-OMEGA+Q/R+R*A+0.5*(1.+Q)* (R*R)*(1.-B*B) +0.5*(1.-Q)
+ DF=A+(1.+Q)*R*(1.-B*B)-Q/(R*R)-(R+A)*SQR3
+ RETURN
+ END
+C
+ SUBROUTINE ROSUB(AX,AY,AZ,BX,BY,BZ,GX,GY,GZ,COSB)
+C ---------------------------------------------------------
+C Used by Roches: makes A and B vectors normal to G,
+C by subtracting their G component, then multiplies by 1/cosb
+C ---------------------------------------------------------
+ COSB1=1./COSB
+ A1=AX*GX+AY*GY+AZ*GZ
+ AX=(AX-GX*A1)*COSB1
+ AY=(AY-GY*A1)*COSB1
+ AZ=(AZ-GZ*A1)*COSB1
+ A1=BX*GX+BY*GY+BZ*GZ
+ BX=(BX-GX*A1)*COSB1
+ BY=(BY-GY*A1)*COSB1
+ BZ=(BZ-GZ*A1)*COSB1
+ RETURN
+ END
+C
+ FUNCTION RTNEWT(FUNCD,X1,X2,XACC,XGUESS,Q,OMEGA,A,B,NERROR)
+C -----------------------------------------------------------
+C Newton-Cotes function's zeroes finder; modified from
+C "Numerical Recipes" By Press,Teukolwsky,Flannery,Vetterling
+C Cambridge University Press 1986
+C -----------------------------------------------------------
+ PARAMETER (JMAX=20)
+C RTNEWT=.5*(X1+X2)
+ RTNEWT=XGUESS
+ DO 11 J=1,JMAX
+ CALL FUNCD(RTNEWT,F,DF,Q,OMEGA,A,B)
+ DX=F/DF
+ RTNEWT=RTNEWT-DX
+ IF((X1-RTNEWT)*(RTNEWT-X2).LT.0.) THEN
+ NERROR=NERROR+1
+ IF(NERROR.LE.10) THEN
+ WRITE(6,1000)
+ 1000 FORMAT(' ERROR ! : Subr. RTNEWT : zero out of bounds!')
+ IF(NERROR.EQ.10) WRITE(6,1100)
+ 1100 FORMAT(' 10 ERRORS !!, ERROR WARNING SUPPRESSED !')
+ ENDIF
+ ENDIF
+ IF(ABS(DX).LT.XACC) RETURN
+ 11 CONTINUE
+ NERROR=NERROR+1
+ IF(NERROR.LE.10) THEN
+ WRITE(6,2000)
+ 2000 FORMAT(' WARNING ! : Subr. RTNEWT : not converging!')
+ IF(NERROR.EQ.10) WRITE(6,1100)
+ ENDIF
+ RETURN
+ END
+C
+ FUNCTION RTSAFE(FUNCD,X1,X2,XACC,Q,DUM,NERROR)
+C ----------------------------------------------------
+C Zeroes for a function FUNCD, from "Numerical Recipes",
+C By Press, Teukolwsky, Flannery, with minor changes
+C ----------------------------------------------------
+ PARAMETER (MAXIT=100)
+ CALL FUNCD(X1,FL,DF,Q,DUM)
+C test if root exactly on bound
+ IF(ABS(FL).LT.XACC) THEN
+ RTSAFE=X1
+ RETURN
+ ENDIF
+ CALL FUNCD(X2,FH,DF,Q,DUM)
+ IF(ABS(FH).LT.XACC) THEN
+ RTSAFE=X2
+ RETURN
+ ENDIF
+ IF(FL*FH.GE.0.) THEN
+ NERROR=NERROR+1
+ IF(NERROR.LE.10) THEN
+ WRITE(6,1000)
+ 1000 FORMAT(' ERROR IN RTSAFE! :root must be bracketed')
+ IF(NERROR.EQ.10) WRITE(6,1100)
+ 1100 FORMAT(' 10 ERRORS !!, ERROR WARNING SUPPRESSED !')
+ ENDIF
+C RTSAFE=0.0
+C RETURN
+ ENDIF
+ IF(FL.LT.0.)THEN
+ XL=X1
+ XH=X2
+ ELSE
+ XH=X1
+ XL=X2
+ SWAP=FL
+ FL=FH
+ FH=SWAP
+ ENDIF
+ RTSAFE=.5*(X1+X2)
+ DXOLD=ABS(X2-X1)
+ DX=DXOLD
+ CALL FUNCD(RTSAFE,F,DF,Q,DUM)
+ DO 11 J=1,MAXIT
+ IF(((RTSAFE-XH)*DF-F)*((RTSAFE-XL)*DF-F).GE.0.
+ * .OR. ABS(2.*F).GT.ABS(DXOLD*DF) ) THEN
+ DXOLD=DX
+ DX=0.5*(XH-XL)
+ RTSAFE=XL+DX
+ IF(XL.EQ.RTSAFE)RETURN
+ ELSE
+ DXOLD=DX
+ DX=F/DF
+ TEMP=RTSAFE
+ RTSAFE=RTSAFE-DX
+ IF(TEMP.EQ.RTSAFE)RETURN
+ ENDIF
+ IF(ABS(DX).LT.XACC) RETURN
+ CALL FUNCD(RTSAFE,F,DF,Q,DUM)
+ IF(F.LT.0.) THEN
+ XL=RTSAFE
+ FL=F
+ ELSE
+ XH=RTSAFE
+ FH=F
+ ENDIF
+ 11 CONTINUE
+ NERROR=NERROR+1
+ IF(NERROR.LE.10) THEN
+ WRITE (6,2000)
+ 2000 FORMAT(' ERROR IN RTSAFE !',
+ 1 ' : Root not found:exceeding maximum iterations')
+ IF(NERROR.EQ.10) WRITE(6,1100)
+ ENDIF
+ RETURN
+ END
+C
+ SUBROUTINE SFERA2(PRINT6,PRINTFILE,NFILE,
+ 1 NTH,NTHF,NCORD,NPUNT,NPUNTC,R,X0,Y0,Z0,
+ 2 NPNTMX,X,Y,Z,G,GX,GY,GZ,FX,FY,FZ,TX,TY,TZ,AREA,A,NUMCOARSE,
+ 3 SINFI,COSFI)
+C ----------------------------------------------------------------
+C Called for a spherical star, fills the surface element list
+C with the surface elements of a sphere of radiur R and center
+C in X0,Y0,Z0. NTH is the number of theta points used to represent
+C the sphere; NFI, the number of phi points, is 2(NTH-1), at equator,
+C approx.: (NTH-1)*SIN(THETA) otherwise, to have equi-area surfaces.
+C Points over and below equator are treated toghether to take
+C advantage of simmetry and shorten computation.
+C
+C COARSE MAP USED FOR REFLECTION COMPUTATION:
+C A coarse map NUMCOARSE is built, containing, for each fine surf.el.
+C the number of its coarse surf. element.
+C To build the coarse map the fine theta points are divided
+C into groups of NTHF points:
+C Each group of NTHF theta values defines a region on the sphere:
+C wide NTHF in theta and 0-2*PI in phi. This region is divided into
+C a number of "PHI" coarse intervals, in such a way to have approx,
+C the same area for coarse surf. elem. of different theta values.
+C The phi values of each theta in this region are divided into
+C this number of coarse intervals.
+C For each theta, in the group of theta values, the coarse phi
+C interval are divided into NTHF fine "PHI" interval, in this way the
+C fine grid surf. element are only approx. equiarea, and each
+C coarse surf. elem. consists of NTHF*NTHF fine elem.
+C
+C ICT is the number of the coarse el. for each fine, it is
+C odd over and even below equator.
+C
+C NFIC is the number of phi coarse for a range of theta
+C
+C NFIF is the number of phi fine points for each coarse, it
+C depends on theta coarse; diffenert theta coarse have a diffenent
+C number of phi.
+C
+C ----------------------------------------------------------------
+ DIMENSION X(NPNTMX),Y(NPNTMX),Z(NPNTMX)
+ DIMENSION G(NPNTMX),GX(NPNTMX),GY(NPNTMX),GZ(NPNTMX)
+ DIMENSION FX(NPNTMX),FY(NPNTMX),FZ(NPNTMX)
+ DIMENSION TX(NPNTMX),TY(NPNTMX),TZ(NPNTMX)
+ DIMENSION A(NPNTMX)
+ DIMENSION NUMCOARSE(NPNTMX)
+ DIMENSION SINFI(NPNTMX),COSFI(NPNTMX)
+ LOGICAL PRINTFILE,PRINT6
+ DATA PI /3.1415926/
+C DATA PI /3.141592653589793/
+C REAL*8 PI,DFI,DTH,FI,TH
+C
+C Delta theta fine: theta has bot end included (fi doesn't)
+ DTH=PI/(NTH-1)
+ DTH5=0.5*DTH
+C Delta phi for equator ( nfi is chosen to have dfi=dth )
+C DFI=(2.*PI)/NFIEQ , DFI5=0.5*DFI
+ NFIEQ=2*(NTH-1)
+C
+C ............... Polar and equatorial interval extension
+C NCORD=1: Theta and phi arch extension for each surface element
+C NCORD=2: Theta and phi cord extension of each surface element
+C NCORD=3: Theta and phi tangent segment of each surface element
+ IF(NCORD.LE.1) THEN
+ RTCORDA=DTH5*R
+ ELSE IF(NCORD.EQ.2) THEN
+ RTCORDA=R*SIN(DTH5)
+ ELSE
+ RTCORDA=R*TAN(DTH5)
+ ENDIF
+C
+C Theta coarse arc interval
+ DTHC=DTH*NTHF
+ DTHC5=DTHC*0.5
+C Number of coarse theta points
+ NTHC=(NTH+(NTHF-1))/NTHF
+ IF(NTHC*NTHF.NE.NTH+NTHF-1) THEN
+ WRITE(6,900) NTH,NTHF
+ IF(PRINTFILE) WRITE(NFILE,900) NTH,NTHF
+ 900 FORMAT(/,' ERRROR! The number of theta points:',I6,
+ 1 ' is not consistent with the coarsing factor:',I6)
+ ENDIF
+C In previous versions NFIF was recomputed at each FI
+C ( Each theta had a different number of phi (fine equiarea)
+ NFIF=NTHF
+C
+ NPUNT=0
+ LOST=0
+ TH=0.0
+C first coarse point number
+ ICT=NPUNTC
+C ---------------
+C ................................ North polar cap
+C ---------------
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 500
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)=0.0
+ Y(NPUNT)=0.0
+ Z(NPUNT)=R
+ GX(NPUNT)=0.0
+ GY(NPUNT)=0.0
+ GZ(NPUNT)=1.0
+ FX(NPUNT)=0.
+ FY(NPUNT)=RTCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)=RTCORDA
+ TY(NPUNT)=0.
+ TZ(NPUNT)=0.
+ A(NPUNT)=(2.*PI)*(1.- COS(DTH5))*R*R
+ NUMCOARSE(NPUNT)=ICT
+C Total area computation
+ AREA=A(NPUNT)
+C ---------------
+C .................................... South polar cap
+C ---------------
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 500
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)=0.0
+ Y(NPUNT)=0.0
+ Z(NPUNT)=-R
+ GX(NPUNT)=0.0
+ GY(NPUNT)=0.0
+ GZ(NPUNT)=-1.0
+ FX(NPUNT)=0.
+ FY(NPUNT)=-RTCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)=RTCORDA
+ TY(NPUNT)=0.
+ TZ(NPUNT)=0.
+ A(NPUNT)=A(NPUNT-1)
+ AREA=AREA+A(NPUNT)
+ NUMCOARSE(NPUNT)=ICT+1
+C
+C ------------------------------
+C ...........................coarse cap: other theta points
+C ------------------------------
+ COSTH2=COS(DTH5)
+ NTHF2=NTHF/2+1
+ DO 5 IT=2,NTHF2
+ TH=DTH*(IT-1)
+ SINTH=SIN(TH)
+ COSTH=COS(TH)
+ COSTH1=COSTH2
+ COSTH2=COS(TH+DTH5)
+C .... Number of phi fines
+C NFI=NINT(NFIEQ*SINTH)
+ NFI=INT(NFIEQ*SINTH)
+ IF(NFI.LT.4) NFI=4
+ DFI=2.*PI/NFI
+ DFI5=DFI*0.5
+ DFIRR=DFI*R*R
+ AREATH=DFIRR*(COSTH1-COSTH2)
+C .... Interval extension for this theta value
+ IF(NCORD.LE.1) THEN
+ RFCORDA=DFI5*R
+C RTCORDA=DTH5*R ! not changing
+ ELSE IF(NCORD.EQ.2) THEN
+ RFCORDA=R*SIN(TH+DTH5)*SIN(DFI5)
+C RTCORDA=R*SIN(DTH5) ! not changing
+ ELSE
+ RFCORDA=R*SIN(TH+DTH5)*TAN(DFI5)/COS(DTH5)
+C RTCORDA=R*TAN(DTH5) ! not changing
+ ENDIF
+C
+ DO 5 IF=1,NFI
+C FI=(IF-1)*DFI
+C For testing purpose i shift phi: ==========================
+C To make phi points consistent with a run: fine=coarse=this nthc
+ FI=DFI * ( IF-1-INT(NFIF/2) )
+ COSFII=COS(FI)
+ SINFII=SIN(FI)
+C
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 5
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ GX(NPUNT)=SINTH*COSFII
+ GY(NPUNT)=SINTH*SINFII
+ GZ(NPUNT)=COSTH
+ X(NPUNT)=R*GX(NPUNT)
+ Y(NPUNT)=R*GY(NPUNT)
+ Z(NPUNT)=R*GZ(NPUNT)
+ FX(NPUNT)=-SINFII*RFCORDA
+ FY(NPUNT)= COSFII*RFCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)= (COSTH*RTCORDA)*COSFII
+ TY(NPUNT)= (COSTH*RTCORDA)*SINFII
+ TZ(NPUNT)=-(SINTH*RTCORDA)
+ A(NPUNT)=AREATH
+ AREA=AREA+AREATH
+ NUMCOARSE(NPUNT)=ICT
+C ........ Simmetric point with theta=pi-theta
+C ( cos(theta) -> - cos(theta) , z -> -z)
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 5
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)= X(NPUNT-1)
+ Y(NPUNT)= Y(NPUNT-1)
+ Z(NPUNT)=-Z(NPUNT-1)
+ GX(NPUNT)= GX(NPUNT-1)
+ GY(NPUNT)= GY(NPUNT-1)
+ GZ(NPUNT)=-GZ(NPUNT-1)
+ FX(NPUNT)=-SINFII*RFCORDA
+ FY(NPUNT)= COSFII*RFCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)=-(COSTH*RTCORDA)*COSFII
+ TY(NPUNT)=-(COSTH*RTCORDA)*SINFII
+ TZ(NPUNT)=-SINTH*RTCORDA
+ A(NPUNT)=AREATH
+ AREA=AREA+AREATH
+ NUMCOARSE(NPUNT)=ICT+1
+ 5 CONTINUE
+C
+C ------------------
+C .................................... other theta points
+C ------------------
+C
+C ------------------------------ Loop on theta values
+C --------------------
+ NTH1=NTHF/2+2
+ NTH2=NTH/2-NTHF/2
+ TH1=(NTH1-1)*DTH
+ NFIC=1
+ COSTH2=COS(TH1-DTH5)
+C
+ DO 20 IT=NTH1,NTH2
+ TH=(IT-1)*DTH
+ SINTH=SIN(TH)
+ COSTH=COS(TH)
+ COSTH1=COSTH2
+ COSTH2=COS(TH+DTH5)
+C
+C ................................
+C ................ CHANGING THE THETA COARSE POINT EACH NTHF FINE POINTS.
+C The first time change when it=nth1, then at each nthf values
+C When NTHF=1 defines all data at each pass. (each IT value)
+C
+ IF(MOD(IT+NTHF/2-1,NTHF).EQ.0) THEN
+C
+C First coarse for next theta group of nthf theta fine points
+C ( by using the preceeding value for ICT)
+ ICT=ICT+NFIC*2
+C
+C Number of phi coarse for this nthf interval
+C THLAST=DTH*(IT-1+NTHF)
+C NFIC=NINT(2.*PI*SIN(THLAST)/DTHC)
+ THMEDIO=DTH*(IT-1+NTHF/2)
+C NFIC=NINT(2.*PI*SIN(THMEDIO)/DTHC)
+ NFIC=INT(2.*PI*SIN(THMEDIO)/DTHC)
+ IF(NFIC.LT.4) NFIC=4
+C
+C Number of phi fine for this coarse theta range
+ NFI=NFIC*NTHF
+ DFI=(2.*PI)/NFI
+ DFI5=DFI*0.5
+ DFIRR=DFI*R*R
+ IF(NCORD.LE.1) THEN
+ RDFI5=R*DFI5
+ ELSE IF(NCORD.EQ.2) THEN
+ RSINDFI5=R*SIN(DFI5)
+ ELSE
+ RTANDFI5=R*TAN(DFI5)
+ ENDIF
+C Cos and sin fi values for this theta range
+ DO 25 IF=1,NFI
+C FI=DFI*(IF-1)
+C For testing purpose: ==========================
+C To make phi points consistent with a run: fine=coarse=this nthc
+ FI=DFI * ( IF-1-INT(NFIF/2) )
+ COSFI(IF)=COS(FI)
+ SINFI(IF)=SIN(FI)
+ 25 CONTINUE
+ ENDIF
+C ................................
+C
+C Set first coarse for the phi loop
+ ICF=ICT
+C Area for the following surface elements
+ AREATH=DFIRR*(COSTH1-COSTH2)
+C ............... interval extension for this theta value
+ IF(NCORD.LE.1) THEN
+ RFCORDA=RDFI5
+ ELSE IF(NCORD.EQ.2) THEN
+ RFCORDA=SIN(TH+DTH5)*RSINDFI5
+C RFCORDA=RSINDFI5
+ ELSE
+ RFCORDA=SIN(TH+DTH5)*RTANDFI5/COS(DTH5)
+C RFCORDA=RTANDFI5
+ ENDIF
+C
+C ------------------------------------ Loop on phi values
+C --------------------
+C
+ DO 20 IF=1,NFI
+C
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 20
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ GX(NPUNT)=SINTH*COSFI(IF)
+ GY(NPUNT)=SINTH*SINFI(IF)
+ GZ(NPUNT)=COSTH
+ X(NPUNT)=R*GX(NPUNT)
+ Y(NPUNT)=R*GY(NPUNT)
+ Z(NPUNT)=R*GZ(NPUNT)
+ FX(NPUNT)=-SINFI(IF)*RFCORDA
+ FY(NPUNT)= COSFI(IF)*RFCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)= (COSTH*RTCORDA)*COSFI(IF)
+ TY(NPUNT)= (COSTH*RTCORDA)*SINFI(IF)
+ TZ(NPUNT)=-(SINTH*RTCORDA)
+ A(NPUNT)=AREATH
+ AREA=AREA+AREATH
+ NUMCOARSE(NPUNT)=ICF
+C ........ Simmetric point with theta=pi-theta
+C ( cos(theta) -> - cos(theta) , z -> -z)
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 20
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)= X(NPUNT-1)
+ Y(NPUNT)= Y(NPUNT-1)
+ Z(NPUNT)=-Z(NPUNT-1)
+ GX(NPUNT)= GX(NPUNT-1)
+ GY(NPUNT)= GY(NPUNT-1)
+ GZ(NPUNT)=-GZ(NPUNT-1)
+ FX(NPUNT)=-SINFI(IF)*RFCORDA
+ FY(NPUNT)= COSFI(IF)*RFCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)=-(COSTH*RTCORDA)*COSFI(IF)
+ TY(NPUNT)=-(COSTH*RTCORDA)*SINFI(IF)
+ TZ(NPUNT)=-SINTH*RTCORDA
+ A(NPUNT)=AREATH
+ AREA=AREA+AREATH
+ NUMCOARSE(NPUNT)=ICF+1
+C
+C AT each NFIF phi points change the coarse phi point
+C Including the remainder of the division IF/NFIF in the last coarse
+ IF(MOD(IF,NFIF).EQ.0) ICF=ICF+2
+C
+ 20 CONTINUE
+C
+C ----------------
+C ................................. Coarse equator
+C ----------------
+ ICT=ICT+NFIC*2
+ NTH1=NTH/2-NTHF/2+1
+ NTH2=NTH/2
+ TH1=DTH*(NTH1-1)
+ COSTH2=COS(TH1-DTH5)
+C Number of phi coarse for equator fine and coarse=NFIEQ
+ NFIC=NINT(2.*PI/DTHC)
+C NFIC=INT(2.*PI/DTHC)
+ IF(NFIC.LT.4) NFIC=4
+ NFI=NFIC*NFIF
+ IF(NFI.NE.NFIEQ) THEN
+ WRITE(6,2000) NFI,NFIEQ
+ IF(PRINTFILE) WRITE(NFILE,2000) NFI,NFIEQ
+ 2000 FORMAT(' ERROR IN SFERA2! NFI,NFIEQ INCONSISTENT!:',I5,1X,I5)
+ ENDIF
+ DFI=(2.*PI)/NFI
+ DFI5=DFI*0.5
+ DFIRR=DFI*R*R
+ DO 27 IF=1,NFI
+C FI=DFI*(IF-1)
+C For testing purpose: =========================
+C To make phi points consistent with a run: fine=coarse=this nthc
+ FI=DFI * ( IF-1-INT(NFIF/2) )
+ COSFI(IF)=COS(FI)
+ SINFI(IF)=SIN(FI)
+ 27 CONTINUE
+C
+ DO 29 IT=NTH1,NTH2
+ TH=(IT-1)*DTH
+ SINTH=SIN(TH)
+ COSTH=COS(TH)
+ COSTH1=COSTH2
+ COSTH2=COS(TH+DTH5)
+C ... Interval extension for this theta value #### QUALCOSA ESTRAIBILE
+ IF(NCORD.LE.1) THEN ! DAL LOOP
+ RFCORDA=DFI5*R
+ ELSE IF(NCORD.EQ.2) THEN
+ RFCORDA=R*SIN(TH+DTH5)*SIN(DFI5)
+C RFCORDA=R*SIN(DFI5)
+ ELSE
+ RFCORDA=R*SIN(TH+DTH5)*TAN(DFI5)/COS(DTH5)
+C RFCORDA=R*TAN(DFI5)
+ ENDIF
+C
+ ICF=ICT
+ AREATH=DFIRR*(COSTH1-COSTH2)
+C ........................ Fi loop
+ DO 29 IF=1,NFI
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 29
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ GX(NPUNT)=SINTH*COSFI(IF)
+ GY(NPUNT)=SINTH*SINFI(IF)
+ GZ(NPUNT)=COSTH
+ X(NPUNT)=R*GX(NPUNT)
+ Y(NPUNT)=R*GY(NPUNT)
+ Z(NPUNT)=R*GZ(NPUNT)
+ FX(NPUNT)=-SINFI(IF)*RFCORDA
+ FY(NPUNT)= COSFI(IF)*RFCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)= (COSTH*RTCORDA)*COSFI(IF)
+ TY(NPUNT)= (COSTH*RTCORDA)*SINFI(IF)
+ TZ(NPUNT)=-(SINTH*RTCORDA)
+ A(NPUNT)=AREATH
+ AREA=AREA+AREATH
+ NUMCOARSE(NPUNT)=ICF
+C ........ Simmetric point with theta=pi-theta
+C ( cos(theta) -> - cos(theta) , z -> -z)
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 29
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ X(NPUNT)= X(NPUNT-1)
+ Y(NPUNT)= Y(NPUNT-1)
+ Z(NPUNT)=-Z(NPUNT-1)
+ GX(NPUNT)= GX(NPUNT-1)
+ GY(NPUNT)= GY(NPUNT-1)
+ GZ(NPUNT)=-GZ(NPUNT-1)
+ FX(NPUNT)=-SINFI(IF)*RFCORDA
+ FY(NPUNT)= COSFI(IF)*RFCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)=-(COSTH*RTCORDA)*COSFI(IF)
+ TY(NPUNT)=-(COSTH*RTCORDA)*SINFI(IF)
+ TZ(NPUNT)=-SINTH*RTCORDA
+ A(NPUNT)=AREATH
+ AREA=AREA+AREATH
+ NUMCOARSE(NPUNT)=ICF
+C ..... At each NFIF phi points change the coarse phi point
+ IF(MOD(IF,NFIF).EQ.0) ICF=ICF+1
+ 29 CONTINUE
+C ------------
+C ......................................... Equator
+C ------------
+ SINTH=1.0
+ COSTH=0.0
+ COSTH1=COS(0.5*PI-DTH5)
+ COSTH2=COS(0.5*PI+DTH5)
+C Cord lenght from last theta value
+ ICF=ICT
+C
+ DO 30 IF=1,NFI
+C
+ IF(NPUNT.EQ.NPNTMX) THEN
+ LOST=LOST+1
+ GOTO 30
+ ELSE
+ NPUNT=NPUNT+1
+ ENDIF
+ GX(NPUNT)=SINTH*COSFI(IF)
+ GY(NPUNT)=SINTH*SINFI(IF)
+ GZ(NPUNT)=COSTH
+ X(NPUNT)=R*GX(NPUNT)
+ Y(NPUNT)=R*GY(NPUNT)
+ Z(NPUNT)=R*GZ(NPUNT)
+ FX(NPUNT)=-SINFI(IF)*RFCORDA
+ FY(NPUNT)= COSFI(IF)*RFCORDA
+ FZ(NPUNT)=0.
+ TX(NPUNT)=COSTH*COSFI(IF)*RTCORDA
+ TY(NPUNT)=COSTH*SINFI(IF)*RTCORDA
+ TZ(NPUNT)=-SINTH*RTCORDA
+ A(NPUNT)=DFI*(COSTH1-COSTH2)*(R*R)
+ AREA=AREA+A(NPUNT)
+ NUMCOARSE(NPUNT)=ICF
+C Change the coarse point at each nfif fine points
+ IF(MOD(IF,NFIF).EQ.0)ICF=ICF+1
+ 30 CONTINUE
+ NPUNTC=ICT+NFIC
+C
+C all complications ended ! WHOW !
+C
+C Shifts star , if requested
+ IF(X0.NE.0.0) THEN
+ DO 50 J=1,NPUNT
+ 50 X(J)=X(J)+X0
+ ENDIF
+ IF(Y0.NE.0.0) THEN
+ DO 53 J=1,NPUNT
+ 53 Y(J)=Y(J)+Y0
+ ENDIF
+ IF(Z0.NE.0.0) THEN
+ DO 56 J=1,NPUNT
+ 56 Z(J)=Z(J)+Z0
+ ENDIF
+C Constant potential
+ DO 59 I=1,NPUNT
+ 59 G(I)=1./R
+C
+ 500 IF(LOST.GT.0) THEN
+ WRITE(6,1000) NPUNT,LOST
+ IF(PRINTFILE) WRITE(NFILE,1000) NPUNT,LOST
+ 1000 FORMAT(' ERROR! OBJECT NOT PROPERLY DRAWN! '/
+ 1 ' space available for only:',I6,' surface points.',
+ 2 ' lost points:',I6/
+ 3 ' REDUCE THE NTHETA PARAMETER OR GO INTO THE FORTRAN LIST TO',
+ 4 ' INCREASE MAXPT')
+ ENDIF
+C
+ RETURN
+ END
+C
+ SUBROUTINE SUMREFL(NPNT,NIF,FBOLREF,FRIFL)
+C --------------------------------------------------
+C This routine sums the reflected surface element flux
+C obtaining the total light reflected by each object
+C ------------------------------------------------
+ DIMENSION FRIFL(NPNT),FBOLREF(3),NIF(2,3)
+C
+C ..................... loop on objects
+ DO 10 I=1,3
+ N1=NIF(1,I)
+ N2=NIF(2,I)
+ FBOLREF(I)=0.0
+ IF(N1.LE.0.OR.N2.LE.0) GOTO 10
+ DO 20 J=N1,N2
+ FBOLREF(I)=FBOLREF(I)+FRIFL(J)
+ 20 CONTINUE
+ 10 CONTINUE
+ RETURN
+ END
+C
+ SUBROUTINE TCHANGE(NPNT,FB,FRIFL,T,A)
+C ----------------------------------------------------
+C Given FB: bolometric flux * area A, computes T,
+C assuming black body spectrum flux=Area*T**4 *ac/4pi
+C ----------------------------------------------------
+ DIMENSION FB(NPNT),T(NPNT),A(NPNT),FRIFL(NPNT)
+ PARAMETER (AC4PI=1.8065E-5)
+ DO 10 I=1,NPNT
+ 10 T(I)=( (FB(I)+FRIFL(I)) / (AC4PI*A(I)) )**0.25
+ RETURN
+ END
+C
+ SUBROUTINE TOTALE(N,T,A)
+C -------------------------------------
+C T=sum of A(i)
+C -------------------------------------
+ DIMENSION A(N)
+ T=0
+ DO 10 I=1,N
+ 10 T=T+A(I)
+ RETURN
+ END
+C
+ SUBROUTINE UPPERC(C,PROMPT,KF)
+C -----------------------------------------------------------
+C Reads an input command and transforms it to upcase letters
+C adding 32 to the ascii character decimal value
+C (32 is the offset between upcase and lowercase letters)
+C This routines doesn't follows the ANSI FORTRAN standards
+C If "?" is encountered the scanning stops and KF=1 is returned
+C IF "$" in first position,the command is a DCL command,
+C and the routine "spawn" creating a subprocess.
+C If "!" in first position the line is skipped being a comment
+C If "!!" in first position this is a title line for printed outp.
+C -----------------------------------------------------------
+ PARAMETER (MAXTITLE=5)
+ COMMON/TITLE/NTITLEMX,NTITLE,TITLE(MAXTITLE)
+ CHARACTER*80 TITLE
+ CHARACTER*(*) C,PROMPT
+C
+ LENGTH=LEN(C)
+C
+ 10 WRITE(6,1000) PROMPT
+ 1000 FORMAT(/A,$)
+C 1000 FORMAT(/' INPUT FLAG > ',$)
+C note:"$": it's a "Vax Fortran Extension" way to give a prompt
+ READ(5,1100,ERR=500,END=600)C
+ 1100 FORMAT(A)
+C ........... vax fortran way to create a subprocess
+ IF(C(1:1).EQ.'$') THEN
+ C=C(2:)
+ ISTATUS=LIB$SPAWN(C)
+ GOTO 10
+C ............ the line is skipped as a comment
+ ELSE IF(C(1:1).EQ.'!') THEN
+C Title line inserted in common
+ IF(C(2:2).EQ.'!'.AND.NTITLE.LT.NTITLEMX) THEN
+ NTITLE=NTITLE+1
+ TITLE(NTITLE)=C(3:)
+ ENDIF
+ GOTO 10
+ ENDIF
+C
+ KF=0
+ DO 20 I=1,LENGTH
+ IF(C(I:I).EQ.'?') THEN
+ KF=1
+ RETURN
+ ENDIF
+ IF(C(I:I).GE.'a'.AND.C(I:I).LE.'z') THEN
+ IC=ICHAR(C(I:I))
+ IC=IC-32
+ C(I:I)=CHAR(IC)
+C 1 C(I:I)=CHAR(ICHAR(C(I:I))+32)
+C 1 C(I:I)=CHAR(IAND( ICHAR(C(I:I)) , '5F'X ) ) Maybe faster
+ ENDIF
+ 20 CONTINUE
+C
+ RETURN
+ 600 WRITE(6,6000)
+ 6000 FORMAT(' ERROR: EOF ENCOUNTERED READING INPUT ')
+ GOTO 10
+ 500 WRITE(6,5000)
+ 5000 FORMAT(' ERROR READING INPUT! Format is A80 Reenter.')
+ GOTO 10
+ END
+C
projects / cbs.git / commitdiff